From a6f23f6a619b913e4114df2f10b90bb8c1a57733 Mon Sep 17 00:00:00 2001
From: Gabriel Gros <gabriel.gros@student-cs.fr>
Date: Tue, 1 Apr 2025 16:17:37 +0200
Subject: [PATCH] =?UTF-8?q?Ajout=20pr=C3=A9conditionnement?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
loss_function.py | 165 +++++++++++++++++++++++++++++++++++++++++++----
1 file changed, 153 insertions(+), 12 deletions(-)
diff --git a/loss_function.py b/loss_function.py
index 59ef63a..40456b7 100644
--- a/loss_function.py
+++ b/loss_function.py
@@ -2,6 +2,7 @@ import numpy as np
import matplotlib.pyplot as plt
import math as mt
from itertools import permutations
+import numdifftools as nd
def generate_data_multi(n, means, scale, weights, distribution='gaussian', kernel=None):
@@ -224,7 +225,7 @@ def batch_algo_gradient_descent(m, n, means, weights, num_steps=100, eta=0.1):
# Dictionary function Phi: for a given tau vector (r,), returns a matrix (N_freq x r)
def Phi(tau):
- return np.exp(2j * np.pi * np.outer(freqs.flatten(), tau))
+ return np.exp(-2j * np.pi * np.outer(freqs.flatten(), tau))
# Model function: Y_model = G * Phi(tau) * a.
def Y_model(a, tau):
@@ -300,7 +301,7 @@ def batch_algo_gradient_descent_without_estimator(m, n, means, weights,est, num_
# Dictionary function Phi: for a given tau vector (r,), returns a matrix (N_freq x r)
def Phi(tau):
- return np.exp(2j * np.pi * np.outer(freqs.flatten(), tau))
+ return np.exp(-2j * np.pi * np.outer(freqs.flatten(), tau))
# Model function: Y_model = G * Phi(tau) * a.
def Y_model(a, tau):
@@ -582,7 +583,7 @@ def show_loss_path(tau_est,tau_star,tau_min,tau_max,nb_values,m,n,num_steps =100
-def show_conv_speed(tau_est,tau_star,m,n,num_steps =1000,eta=0.1):
+def show_conv_speed_tau(tau_est,tau_star,m,n,num_steps =1000,eta=0.1):
"""
We take r=2, and let a1 = a2 = 1/2.
@@ -670,8 +671,6 @@ def show_conv_speed(tau_est,tau_star,m,n,num_steps =1000,eta=0.1):
grad_tau = np.real(( (a_est.conj()) * (( Lambda @ G_mat @ Phi(tau_est)).conj().T) ) @ r_vec)
-
-
# Preconditioner for tau-update:
P_tau = np.diag(1/(-Kpp0) * ((a_est.flatten())**2))
#P_tau = 1
@@ -706,16 +705,15 @@ def show_conv_speed(tau_est,tau_star,m,n,num_steps =1000,eta=0.1):
-def precond_diag(a_k,freqs,r,n):
+def precond_diag(a_k,tau_k,freqs,r,n):
"""
Build the adaptative diagonal preconditionner matrix.
"""
K_second_0 = 0
- N = 2*n+1
- L=fourier(freqs).flatten()
- for k in range(N):
- K_second_0 += (abs(L[k])**2) * (k-n)/N
- K_second_0 = (-4*np.pi**2)*K_second_0
+
+ for f in freqs.flatten():
+ K_second_0 += (abs(fourier(f))**2) * (f**2) * (-4*(np.pi**2))
+
P = np.eye(2*r)
@@ -725,9 +723,152 @@ def precond_diag(a_k,freqs,r,n):
return P
-def precond_full(a_k,freqs,r,n):
+def precond_full(a_k,tau_k,freqs,r,n):
"""
Build the adaptative full preconditionner matrix.
"""
+ N=2*n+1
+
Ak = np.diag(np.concatenate((np.array([1 for _ in range(r)]), np.array(a_k))))
+
+
+ G_diag = fourier(freqs).flatten() # shape (N_freq,)
+ G_mat = np.diag(G_diag) # shape (N_freq, N_freq)
+
+ Phi= np.exp(-2j * np.pi * np.outer(freqs.flatten(), tau_k))
+
+
+
+ diag_elements = -2j * np.pi * (np.arange(N) - n) / N
+ Lambda = np.diag(diag_elements)
+
+
+ W = np.hstack((G_mat @ Phi , Lambda @ G_mat @ Phi))
+
+ P = np.linalg.inv( Ak.conj().T @ ((W.conj().T) @ W) @ Ak)
+
+ return P
+
+
+def matrix_S(freqs,a_star,n,r):
+ """
+ Buils the matrix S that adimensionnate the space : With S, the convergence of the gradient descent is linear.
+ """
+ K_second_0 = 0
+
+ for f in freqs.flatten():
+ K_second_0 += (abs(fourier(f))**2) * (f**2) * (-4*(np.pi**2))
+
+ S = np.diag(np.concatenate(([1/k for k in a_star] , [np.sqrt(-K_second_0) for _ in range(r)])))
+
+ return S
+
+
+
+
+def show_conv_speed(theta_est,theta_star,n,precond,num_steps =2000,eta=0.1):
+ """
+ We fix a preconditionner precond, and we look at the convergence of the gradient descent wrt to this preconditionner.
+ We assume here the data is perfect : m is useless
+ """
+ assert len(theta_star) == len(theta_est) , "Error: star and est not same size"
+ assert len(theta_star)%2 ==0 , "Error : r is not even"
+
+ r = len(theta_star)//2
+ freqs=sample_uniform_frequencies(n)
+
+ #Make the input np.array
+ theta_star = np.array(theta_star)
+ theta_est = np.array(theta_est)
+
+ tau_star = theta_star[r:]
+ a_star = theta_star[:r]
+ N=2*n+1
+
+ tau_est = theta_est[r:]
+ a_est = theta_est[:r]
+
+
+ #Reshape the a vectors to make computations easier
+ a_est = a_est.reshape(-1, 1)
+ a_star = a_star.reshape(-1,1)
+
+ S = matrix_S(freqs,a_star.flatten(),n,r)
+ #vitesse_conv = [permute_and_compute_norm(theta_est,theta_star,r)[0]]
+
+ vitesse_conv = [permute_and_compute_norm( (S @ theta_est.reshape(-1,1)).flatten() ,(S @ theta_est.reshape(-1,1)).flatten() ,r)[0]]
+
+ # Precompute the "G" matrix: diagonal with entries given by kernel's Fourier transform at freqs.
+ G_diag = fourier(freqs).flatten() # shape (N_freq,)
+ G_mat = np.diag(G_diag) # shape (N_freq, N_freq)
+
+ # Dictionary function Phi: for a given tau vector (r,), returns a matrix (N_freq x r)
+ def Phi(tau):
+ return np.exp(-2j * np.pi * np.outer(freqs.flatten(), tau))
+
+ # Model function: Y_model = G * Phi(tau) * a.
+ def Y_model(a, tau):
+ return G_mat @ (Phi(tau) @ a)
+
+ Y = Y_model(a_star,tau_star) # No noise assumed, perfect data m -> +inf
+
+ # Run batch gradient descent updates
+ for step in range(num_steps):
+ # Compute the current model prediction
+ Y_estimated = Y_model(a_est, tau_est)
+
+ # Residual: difference between model prediction and observed FCE
+ r_vec = Y_estimated - Y # shape (N_freq, 1)
+
+ # Gradient with respect to amplitudes: grad_a = (G*Phi(tau))^H * r
+ grad_a = (G_mat @ Phi(tau_est)).conj().T @ r_vec # shape (r, 1)
+
+ # Gradient with respect to tau:
+ diag_elements = -2j * np.pi * (np.arange(N) - n) / N
+ Lambda = np.diag(diag_elements)
+
+ grad_tau = np.real(( (a_est.conj()) * (( Lambda @ G_mat @ Phi(tau_est)).conj().T) ) @ r_vec)
+
+ #Concatenation of both gradients
+ grad = np.vstack((grad_a,grad_tau))
+
+ def Loss_theta(theta):
+ return loss(theta,theta_star,r,freqs)
+
+ # Gradient numérique si besoin
+ #grad2 = nd.Gradient(Loss_theta)(theta_est)
+
+ #Preconditionner for the update
+ P = precond(a_est.flatten(),tau_est,freqs,r,n)
+
+ # Gradient descent update with learning rate eta
+ theta_est = theta_est - eta * (P @ grad).flatten()
+ ##theta_est = theta_est -eta * grad.flatten()
+ a_est = theta_est[:r].reshape(-1,1)
+ tau_est = theta_est[r:]
+
+
+ theta_est = np.concatenate((a_est.flatten(),tau_est))
+
+
+ #vitesse_conv.append(permute_and_compute_norm(theta_est,theta_star,r)[0])
+ vitesse_conv.append(permute_and_compute_norm( (S @ theta_est.reshape(-1,1)).flatten() ,(S @ theta_est.reshape(-1,1)).flatten() ,r)[0])
+
+
+ # Final estimated parameters (concatenated amplitudes and centroids)
+ theta_est_final = np.concatenate((a_est.flatten(), tau_est))
+
+ norm , theta_est_final_prime = permute_and_compute_norm(theta_est_final,theta_star,r)
+
+
+ print(vitesse_conv)
+ plt.plot( vitesse_conv , label = f"{precond}")
+ plt.yscale('log')
+ plt.xlabel("Itération")
+ plt.ylabel(f"Norme de la fonction de perte (échelle log)")
+ plt.legend()
+ plt.show()
+
+ return norm, theta_est_final , theta_est_final_prime
+
--
GitLab