矩阵方程求解

稀疏矩阵方程求解：优化算法

正交匹配追踪法

贪婪算法 不求整体最优解，而是试图尽快找到在某种意义上的局部最优解。
典型的贪婪算法有以下匹配追踪算法：
（1）匹配追踪（matching pursuit,MP）法
  基本思想是，不是针对某个代价函数进行最小化，而是考虑迭代地构造一个稀疏解$x:$只使用字典矩阵$\Phi$的少数列向量的线性组合对观测向量$x$实现稀疏逼近$\Phi x=y$，其中字典矩阵$\Phi$被选择的列向量所组成的作用集是以逐列的的方式建立的。在每一步迭代，字典矩阵中同当前残差向量$r=\Phi x-y$最相似的列向量被选择作为作用集的新的一列。如果残差随着迭代的进行递减，则可以保证算法收敛。
代码如下：

function [x_rec] = MatchingPursuit( y, K, Theta)
%% 匹配追踪算法
%% Vsesion:1.0 Written by zhenhuaLiu@ 2021.11.16 HIT ATCI
%% 输入参数：
    %Theta:感知矩阵
    %y: 压缩采样数据
    %K: 稀疏度
%% 输出参数：
    %x_rex: 重构的x向量
A = Theta;
[m, n] = size(Theta);
r = y;%初始化残差
x_rec = zeros(n, 1);%重构的稀疏向量初始化为0向量 
A = Theta;
position = zeros(K, 1);%最大内积的索引值
max_inner = zeros(K, 1);%最大的内积值
for i = 1:K
   % i = 2;
    inner_product = r'*A;%计算残差与感知矩阵的内积；
    [~, max_index] = max(abs(inner_product));% 找到内积的绝对值的最大值的下标
    max_inner(i) = inner_product(max_index);%获得最大的内积值
    r = r - max_inner(i)*A(:,max_index);
    A(:, max_index) = zeros(m,1);%将支撑集索引对用的感知矩阵的该列全部置零
    position(i) = max_index;
end   
x_rec(position) = max_inner;  %重构的稀疏向量除了索引值位置之外其他位置都为0

（2）正交匹配追踪（orthogonal matching pursuit,OMP）法
  匹配追踪只能保证残差向量与每一步迭代所选择的字典矩阵列向量正交，但与以前选择的列向量一般不正交。正交匹配追踪则能够保证每步迭代后残差向量与以前选择的所有列向量正交，以保证迭代的最优性，从而减少了迭代次数，性能也更稳健。正交匹配追踪算法复杂度为O(mn)，可以得到稀疏度$K⩽m/(2\log n)$的系数向量。

import numpy as np

def omp(y, D, S=20):
    M, N=D.shape[0],D.shape[1]
    result = np.zeros((N, 1))
    indices = []
    R = y
    rr=0
    for i in range(S):
        projection = np.dot(D.T, R)
        pos = projection.index(max(projection))
        indices.append(pos)
        At = D[:,indices[:]]
        Ls = np.dot(np.linalg.pinv(At),y)
        R = y - np.dot(At, Ls)
        rr=R
        if (R**2).sum()<1e-6:
            break
    for t, s in zip(indices, Ls):
        result[t] = s
    return result, rr

（3）正则正交匹配追踪（ROMP）法
  在OMP算法基础上，加入正则化过程。首先根据相关原子挑选多个原子作为候选集，然后从候选集中按照正则化原则挑选出部分原子，最后将其并入最终的支撑集，实现原子的快速、有效选择。
代码：

function [ theta ] = CS_ROMP( y,A,K )
%CS_ROMP Summary of this function goes here
%Version: 1.0 written by jbb0523 @2015-04-24
%   Detailed explanation goes here
%   y = Phi * x
%   x = Psi * theta
%   y = Phi*Psi * theta
%   令 A = Phi*Psi, 则y=A*theta
%   现在已知y和A，求theta
%   Reference:Needell D，Vershynin R．Signal recovery from incomplete and
%   inaccurate measurements via regularized orthogonal matching pursuit[J]．
%   IEEE Journal on Selected Topics in Signal Processing，2010，4(2)：310—316.
    [y_rows,y_columns] = size(y);
    if y_rows<y_columns
        y = y';%y should be a column vector
    end
    [M,N] = size(A);%传感矩阵A为M*N矩阵
    theta = zeros(N,1);%用来存储恢复的theta(列向量)
    At = zeros(M,3*K);%用来迭代过程中存储A被选择的列
    Pos_theta = zeros(1,2*K);%用来迭代过程中存储A被选择的列序号
    Index = 0;
    r_n = y;%初始化残差(residual)为y
    %Repeat the following steps K times(or until |I|>=2K)
    for ii=1:K%迭代K次
        product = A'*r_n;%传感矩阵A各列与残差的内积
        %[val,pos] = max(abs(product));%找到最大内积绝对值，即与残差最相关的列
        [val,pos] = Regularize(product,K);%按正则化规则选择原子
        At(:,Index+1:Index+length(pos)) = A(:,pos);%存储这几列
        Pos_theta(Index+1:Index+length(pos)) = pos;%存储这几列的序号
        if Index+length(pos)<=M%At的行数大于列数，此为最小二乘的基础(列线性无关)
            Index = Index+length(pos);%更新Index，为下次循环做准备
        else%At的列数大于行数，列必为线性相关的,At(:,1:Index)'*At(:,1:Index)将不可逆
            break;%跳出for循环
        end
        A(:,pos) = zeros(M,length(pos));%清零A的这几列(其实此行可以不要因为它们与残差正交)
        %y=At(:,1:Index)*theta，以下求theta的最小二乘解(Least Square)
        theta_ls = (At(:,1:Index)'*At(:,1:Index))^(-1)*At(:,1:Index)'*y;%最小二乘解
        %At(:,1:Index)*theta_ls是y在At(:,1:Index)列空间上的正交投影
        r_n = y - At(:,1:Index)*theta_ls;%更新残差
        if norm(r_n)<1e-6%Repeat the steps until r=0
            break;%跳出for循环
        end
        if Index>=2*K%or until |I|>=2K
            break;%跳出for循环
        end
    end
    theta(Pos_theta(1:Index))=theta_ls;%恢复出的theta
end