Fondements de l'algorithme GPSR

L'algorithme de Projection de Gradient pour la Reconstruction Sparse (GPSR) constitue une méthode robuste pour résoudre les problèmes d'optimisation convexe rencontrés en acquisition comprimée (Compressed Sensing). L'objectif principal est de récupérer un signal parsimonieux à partir d'un nombre réduit de mesures linéaires bruitées.

Le modèle mathématique standard s'appuie sur la minimisation de la norme L1, formulée comme suit :

min 0.5 * ||Ax - y||₂² + λ||x||₁

Où A représente la matrice de mesure, y le vecteur d'observation, et λ le paramètre de régularisation qui équilibre la fidélité des données et la parcimonie de la solution.

Implémentation de la fonction GPSR

Contrairement aux méthodes de programmation linéaire classiques, le GPSR transforme le problème initial en un problème de programmation quadratique avec contraintes de boîte ou utilise des techniques de séparation de variables pour accélérer la convergence.

function [sol, stats] = gpsr_solver(Phi, b, tau, limit_iter, stop_crit)
    % Phi : Matrice de détection (M x N)
    % b   : Vecteur de mesures (M x 1)
    % tau : Coefficient de régularisation
    
    [~, N] = size(Phi);
    x_current = zeros(N, 1);
    stats.history = [];

    % Pré-calcul pour optimisation
    PhiT_b = Phi' * b;

    for k = 1:limit_iter
        % Calcul du résidu et du gradient
        residual = Phi * x_current - b;
        grad = Phi' * residual;
        
        % Mise à jour par seuillage doux (Proximal Operator)
        x_prev = x_current;
        alpha = 0.01; % Pas d'apprentissage fixe ou adaptatif
        
        % Descente de gradient projetée
        temp_x = x_prev - alpha * grad;
        x_current = sign(temp_x) .* max(abs(temp_x) - tau * alpha, 0);
        
        % Calcul de la fonction objectif
        current_obj = 0.5 * norm(Phi * x_current - b)^2 + tau * norm(x_current, 1);
        stats.history = [stats.history, current_obj];
        
        % Vérification de la convergence
        if norm(x_current - x_prev) / norm(x_prev + eps) < stop_crit
            fprintf('Convergence atteinte à l''itération %d\n', k);
            break;
        end
    end
    sol = x_current;
end

Simulation et Validation

Pour valider l'efficacité de la reconstruction, nous générons un signal synthétique parsimonieux dans le domaine fréquentiel (DCT) et appliquons une matrice de mesure Gaussienne.

% Configuration de l'expérience
n_size = 1024;       % Dimension du signal
s_level = 40;        % Nombre d'éléments non nuls
m_measure = 300;     % Nombre de mesures

% Création du signal source
original_signal = zeros(n_size, 1);
indices = randperm(n_size, s_level);
original_signal(indices) = randn(s_level, 1);

% Matrice de projection et acquisition
ProjectionMat = randn(m_measure, n_size) / sqrt(m_measure);
observations = ProjectionMat * original_signal + 0.01 * randn(m_measure, 1);

% Exécution de la reconstruction
reg_val = 0.1;
[reconstructed_x, log_data] = gpsr_solver(ProjectionMat, observations, reg_val, 1000, 1e-6);

% Évaluation des performances (SNR)
signal_noise_ratio = 20 * log10(norm(original_signal) / norm(original_signal - reconstructed_x));
fprintf('Rapport Signal sur Bruit (SNR) : %.2f dB\n', signal_noise_ratio);

Extensions Algorithmiques

L'approche GPSR peut être étendue pour répondre à des besoins spécifiques de performance ou de précision :

• Régularisation Lp (0 < p < 1) : Pour induire une parcimonie plus agressive, bien que le problème devienne non-convexe.
• Accélération GPU : L'utilisation de types gpuArray dans MATLAB permet de déporter les multiplications matricielles Phi * x sur le processeur graphique pour traiter des données massives.
• Variantes BB (Barzilai-Borwein) : Utilisation d'un pas adaptatif basé sur la courbure locale pour réduire drastiquement le nombre d'itérations.

% Exemple rapide d'accélération GPU
if canUseGPU
    Phi_gpu = gpuArray(single(ProjectionMat));
    b_gpu = gpuArray(single(observations));
    [sol_gpu, ~] = gpsr_solver(Phi_gpu, b_gpu, reg_val, 500, 1e-5);
    final_result = gather(sol_gpu);
end

Le choix du paramètre λ reste critique : une valeur trop élevée risque de supprimer des composantes essentielles du signal, tandis qu'une valeur trop faible laisse subsister un bruit important dans la solution reconstruite.