myhomotopy.m

来自「code for homotopy in compressed sensing」· M 代码 · 共 50 行

function [x, iterIdx] = myHomotopy(A,y)

iterTimes = 1000;

n = size(A, 1);
m = size(A, 2);
x = zeros(m,1);
actSet = [];

for iterIdx = 1 : iterTimes
    % compute residual correlations
    c = A'*(y-A*x);
    
    % compute active set 
    [lambda, maxIdx] = max(abs(c));   
    actSet = find(abs(abs(c) - lambda) < 1e-5);
    state  = zeros(1,m);    
    state(actSet) = 1;
    
    % compute direction
    R = A(:, actSet)'*A(:,actSet);
    d = inv(R)*sign(c(actSet));
    
    % compute step
    gamma = 1000;
    for idx = 1 : m
        if state(idx)
            % active elements
            myId = find(actSet==idx);
            tmp = max(0,-x(idx)/d(myId));
        else
            % null elements
            av = A(:,idx)'*A(:,actSet)*d;            
            tmp = max(0,(lambda-c(idx)) / (1-av));
            tmp = min(tmp, max(0,(lambda+c(idx)) / (1+av)));
        end
        
        if tmp > 0
           gamma = min(tmp, gamma);
        end
    end
    
    % jump to the next point
    x(actSet) = x(actSet) + gamma*d;
    
    % is it done?
    if norm(y-A*x) < 1e-6
        break;
    end
end