Mercurial > hg > pycsalgos

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/matlab/GAP/GAP.m	Sun Nov 06 20:58:11 2011 +0000
@@ -0,0 +1,144 @@
+function [xhat, Lambdahat] = GAP(y, M, MH, Omega, OmegaH, params, xinit)
+
+%%
+% [xhat, Lambdahat] = GAP(y, M, MH, Omega, OmegaH, params, xinit)
+%
+% Greedy Analysis Pursuit Algorithm
+% This aims to find an approximate (sometimes exact) solution of
+%    xhat = argmin || Omega * x ||_0   subject to   || y - M * x ||_2 <= epsilon.
+%
+% Outputs:
+%   xhat : estimate of the target cosparse vector x0.
+%   Lambdahat : estimate of the cosupport of x0.
+%
+% Inputs:
+%   y : observation/measurement vector of a target cosparse solution x0,
+%       given by relation  y = M * x0 + noise.
+%   M : measurement matrix. This should be given either as a matrix or as a function handle
+%       which implements linear transformation.
+%   MH : conjugate transpose of M.
+%   Omega : analysis operator. Like M, this should be given either as a matrix or as a function
+%           handle which implements linear transformation.
+%   OmegaH : conjugate transpose of OmegaH.
+%   params : parameters that govern the behavior of the algorithm (mostly).
+%      params.num_iteration : GAP performs this number of iterations.
+%      params.greedy_level : determines how many rows of Omega GAP eliminates at each iteration.
+%                            if the value is < 1, then the rows to be eliminated are determined by
+%                                j : |omega_j * xhat| > greedy_level * max_i |omega_i * xhat|.
+%                            if the value is >= 1, then greedy_level is the number of rows to be
+%                            eliminated at each iteration.
+%      params.stopping_coefficient_size : when the maximum analysis coefficient is smaller than
+%                                         this, GAP terminates.
+%      params.l2solver : legitimate values are 'pseudoinverse' or 'cg'. determines which method
+%                        is used to compute
+%                        argmin || Omega_Lambdahat * x ||_2   subject to  || y - M * x ||_2 <= epsilon.
+%      params.l2_accuracy : when l2solver is 'cg', this determines how accurately the above
+%                           problem is solved.
+%      params.noise_level : this corresponds to epsilon above.
+%   xinit : initial estimate of x0 that GAP will start with. can be zeros(d, 1).
+%
+% Examples:
+%
+% Not particularly interesting:
+% >> d = 100; p = 110; m = 60;
+% >> M = randn(m, d);
+% >> Omega = randn(p, d);
+% >> y = M * x0 + noise;
+% >> params.num_iteration = 40;
+% >> params.greedy_level = 0.9;
+% >> params.stopping_coefficient_size = 1e-4;
+% >> params.l2solver = 'pseudoinverse';
+% >> [xhat, Lambdahat] = GAP(y, M, M', Omega, Omega', params, zeros(d, 1));
+%
+% Assuming that FourierSampling.m, FourierSamplingH.m, FDAnalysis.m, etc. exist:
+% >> n = 128;
+% >> M = @(t) FourierSampling(t, n);
+% >> MH = @(u) FourierSamplingH(u, n);
+% >> Omega = @(t) FDAnalysis(t, n);
+% >> OmegaH = @(u) FDSynthesis(t, n);
+% >> params.num_iteration = 1000;
+% >> params.greedy_level = 50;
+% >> params.stopping_coefficient_size = 1e-5;
+% >> params.l2solver = 'cg';   % in fact, 'pseudoinverse' does not even make sense.
+% >> [xhat, Lambdahat] = GAP(y, M, MH, Omega, OmegaH, params, zeros(d, 1));
+%
+% Above: FourierSampling and FourierSamplingH are conjugate transpose of each other.
+%        FDAnalysis and FDSynthesis are conjugate transpose of each other.
+%        These routines are problem specific and need to be implemented by the user.
+
+d = length(xinit(:));
+
+if strcmp(class(Omega), 'function_handle')
+    p = length(Omega(zeros(d,1)));
+else    %% Omega is a matrix
+    p = size(Omega, 1);
+end
+
+iter = 0;
+lagmult = 1e-4;
+Lambdahat = 1:p;
+while iter < params.num_iteration
+    iter = iter + 1;
+    [xhat, analysis_repr, lagmult] = ArgminOperL2Constrained(y, M, MH, Omega, OmegaH, Lambdahat, xinit, lagmult, params);
+    [to_be_removed, maxcoef] = FindRowsToRemove(analysis_repr, params.greedy_level);
+    %disp(['** maxcoef=', num2str(maxcoef), ' target=', num2str(params.stopping_coefficient_size), ' rows_remaining=', num2str(length(Lambdahat)), ' lagmult=', num2str(lagmult)]);
+    if check_stopping_criteria(xhat, xinit, maxcoef, lagmult, Lambdahat, params)
+        break;
+    end
+    xinit = xhat;
+    Lambdahat(to_be_removed) = [];
+
+    %n = sqrt(d);
+    %figure(9);
+    %RR = zeros(2*n, n-1);
+    %RR(Lambdahat) = 1;
+    %XD = ones(n, n);
+    %XD(:, 2:end) = XD(:, 2:end) .* RR(1:n, :);
+    %XD(:, 1:(end-1)) = XD(:, 1:(end-1)) .* RR(1:n, :);
+    %XD(2:end, :) = XD(2:end, :) .* RR((n+1):(2*n), :)';
+    %XD(1:(end-1), :) = XD(1:(end-1), :) .* RR((n+1):(2*n), :)';
+    %XD = FD2DiagnosisPlot(n, Lambdahat);
+    %imshow(XD);
+    %figure(10);
+    %imshow(reshape(real(xhat), n, n));
+end
+return;
+
+
+function [to_be_removed, maxcoef] = FindRowsToRemove(analysis_repr, greedy_level)
+
+    abscoef = abs(analysis_repr(:));
+    n = length(abscoef);
+    maxcoef = max(abscoef);
+    if greedy_level >= 1
+        qq = quantile(abscoef, 1-greedy_level/n);
+    else
+        qq = maxcoef*greedy_level;
+    end
+
+    to_be_removed = find(abscoef >= qq);
+    return;
+
+function r = check_stopping_criteria(xhat, xinit, maxcoef, lagmult, Lambdahat, params)
+
+    r = 0;
+
+    if isfield(params, 'stopping_coefficient_size') && maxcoef < params.stopping_coefficient_size
+        r = 1;
+        return;
+    end
+
+    if isfield(params, 'stopping_lagrange_multiplier_size') && lagmult > params.stopping_lagrange_multiplier_size
+        r = 1;
+        return;
+    end
+
+    if isfield(params, 'stopping_relative_solution_change') && norm(xhat-xinit)/norm(xhat) < params.stopping_relative_solution_change
+        r = 1;
+        return;
+    end
+
+    if isfield(params, 'stopping_cosparsity') && length(Lambdahat) < params.stopping_cosparsity
+        r = 1;
+        return;
+    end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pyCSalgos/GAP/gap.py	Sun Nov 06 20:58:11 2011 +0000
@@ -0,0 +1,312 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Oct 13 14:05:22 2011
+
+@author: ncleju
+"""
+
+#from numpy import *
+#from scipy import *
+import numpy as np
+import scipy as sp
+from scipy import linalg
+import math
+
+from numpy.random import RandomState
+rng = RandomState()
+
+
+
+def Generate_Analysis_Operator(d, p):
+  # generate random tight frame with equal column norms
+  if p == d:
+    T = rng.randn(d,d);
+    [Omega, discard] = np.qr(T);
+  else:
+    Omega = rng.randn(p, d);
+    T = np.zeros((p, d));
+    tol = 1e-8;
+    max_j = 200;
+    j = 1;
+    while (sum(sum(abs(T-Omega))) > np.dot(tol,np.dot(p,d)) and j < max_j):
+        j = j + 1;
+        T = Omega;
+        [U, S, Vh] = sp.linalg.svd(Omega);
+        V = Vh.T
+        #Omega = U * [eye(d); zeros(p-d,d)] * V';
+        Omega2 = np.dot(np.dot(U, np.concatenate((np.eye(d), np.zeros((p-d,d))))), V.transpose())
+        #Omega = diag(1./sqrt(diag(Omega*Omega')))*Omega;
+        Omega = np.dot(np.diag(1 / np.sqrt(np.diag(np.dot(Omega2,Omega2.transpose())))), Omega2)
+    #end
+    ##disp(j);
+#end
+  return Omega
+
+
+def Generate_Data_Known_Omega(Omega, d,p,m,k,noiselevel, numvectors, normstr):
+  #function [x0,y,M,LambdaMat] = Generate_Data_Known_Omega(Omega, d,p,m,k,noiselevel, numvectors, normstr)
+
+  # Building an analysis problem, which includes the ingredients:
+  #   - Omega - the analysis operator of size p*d
+  #   - M is anunderdetermined measurement matrix of size m*d (m<d)
+  #   - x0 is a vector of length d that satisfies ||Omega*x0||=p-k
+  #   - Lambda is the true location of these k zeros in Omega*x0
+  #   - a measurement vector y0=Mx0 is computed
+  #   - noise contaminated measurement vector y is obtained by
+  #     y = y0 + n where n is an additive gaussian noise with norm(n,2)/norm(y0,2) = noiselevel
+  # Added by Nic:
+  #   - Omega = analysis operator
+  #   - normstr: if 'l0', generate l0 sparse vector (unchanged). If 'l1',
+  #   generate a vector of Laplacian random variables (gamma) and
+  #   pseudoinvert to find x
+
+  # Omega is known as input parameter
+  #Omega=Generate_Analysis_Operator(d, p);
+  # Omega = randn(p,d);
+  # for i = 1:size(Omega,1)
+  #     Omega(i,:) = Omega(i,:) / norm(Omega(i,:));
+  # end
+
+  #Init
+  LambdaMat = np.zeros((k,numvectors))
+  x0 = np.zeros((d,numvectors))
+  y = np.zeros((m,numvectors))
+
+  M = rng.randn(m,d);
+
+  #for i=1:numvectors
+  for i in range(0,numvectors):
+
+    # Generate signals
+    #if strcmp(normstr,'l0')
+    if normstr == 'l0':
+        # Unchanged
+
+        #Lambda=randperm(p);
+        Lambda = rng.permutation(int(p));
+        Lambda = np.sort(Lambda[0:k]);
+        LambdaMat[:,i] = Lambda; # store for output
+
+        # The signal is drawn at random from the null-space defined by the rows
+        # of the matreix Omega(Lambda,:)
+        [U,D,Vh] = sp.linalg.svd(Omega[Lambda,:]);
+        V = Vh.T
+        NullSpace = V[:,k:];
+        #print np.dot(NullSpace, rng.randn(d-k,1)).shape
+        #print x0[:,i].shape
+        x0[:,i] = np.squeeze(np.dot(NullSpace, rng.randn(d-k,1)));
+        # Nic: add orthogonality noise
+        #     orthonoiseSNRdb = 6;
+        #     n = randn(p,1);
+        #     #x0(:,i) = x0(:,i) + n / norm(n)^2 * norm(x0(:,i))^2 / 10^(orthonoiseSNRdb/10);
+        #     n = n / norm(n)^2 * norm(Omega * x0(:,i))^2 / 10^(orthonoiseSNRdb/10);
+        #     x0(:,i) = pinv(Omega) * (Omega * x0(:,i) + n);
+
+    #elseif strcmp(normstr, 'l1')
+    elif normstr == 'l1':
+        print('Nic says: not implemented yet')
+        raise Exception('Nic says: not implemented yet')
+        #gamma = laprnd(p,1,0,1);
+        #x0(:,i) = Omega \ gamma;
+    else:
+        #error('normstr must be l0 or l1!');
+        print('Nic says: not implemented yet')
+        raise Exception('Nic says: not implemented yet')
+    #end
+
+    # Acquire measurements
+    y[:,i]  = np.dot(M, x0[:,i])
+
+    # Add noise
+    t_norm = np.linalg.norm(y[:,i],2);
+    n = np.squeeze(rng.randn(m, 1));
+    y[:,i] = y[:,i] + noiselevel * t_norm * n / np.linalg.norm(n, 2);
+  #end
+
+  return x0,y,M,LambdaMat
+
+#####################
+
+#function [xhat, arepr, lagmult] = ArgminOperL2Constrained(y, M, MH, Omega, OmegaH, Lambdahat, xinit, ilagmult, params)
+def ArgminOperL2Constrained(y, M, MH, Omega, OmegaH, Lambdahat, xinit, ilagmult, params):
+
+    #
+    # This function aims to compute
+    #    xhat = argmin || Omega(Lambdahat, :) * x ||_2   subject to  || y - M*x ||_2 <= epsilon.
+    # arepr is the analysis representation corresponding to Lambdahat, i.e.,
+    #    arepr = Omega(Lambdahat, :) * xhat.
+    # The function also returns the lagrange multiplier in the process used to compute xhat.
+    #
+    # Inputs:
+    #    y : observation/measurements of an unknown vector x0. It is equal to M*x0 + noise.
+    #    M : Measurement matrix
+    #    MH : M', the conjugate transpose of M
+    #    Omega : analysis operator
+        #    OmegaH : Omega', the conjugate transpose of Omega. Also, synthesis operator.
+    #    Lambdahat : an index set indicating some rows of Omega.
+    #    xinit : initial estimate that will be used for the conjugate gradient algorithm.
+    #    ilagmult : initial lagrange multiplier to be used in
+    #    params : parameters
+    #        params.noise_level : this corresponds to epsilon above.
+    #        params.max_inner_iteration : `maximum' number of iterations in conjugate gradient method.
+    #        params.l2_accurary : the l2 accuracy parameter used in conjugate gradient method
+    #        params.l2solver : if the value is 'pseudoinverse', then direct matrix computation (not conjugate gradient method) is used. Otherwise, conjugate gradient method is used.
+    #
+
+    #d = length(xinit)
+    d = xinit.size
+    lagmultmax = 1e5;
+    lagmultmin = 1e-4;
+    lagmultfactor = 2;
+    accuracy_adjustment_exponent = 4/5;
+    lagmult = max(min(ilagmult, lagmultmax), lagmultmin);
+    was_infeasible = 0;
+    was_feasible = 0;
+
+    #######################################################################
+    ## Computation done using direct matrix computation from matlab. (no conjugate gradient method.)
+    #######################################################################
+    #if strcmp(params.l2solver, 'pseudoinverse')
+    if params['solver'] == 'pseudoinverse':
+    #if strcmp(class(M), 'double') && strcmp(class(Omega), 'double')
+      if M.dtype == 'float64' and Omega.dtype == 'double':
+        while 1:
+            alpha = math.sqrt(lagmult);
+            #xhat = np.concatenate((M, alpha*Omega(Lambdahat,:)]\[y; zeros(length(Lambdahat), 1)];
+            xhat = np.concatenate((M, np.linalg.lstsq(alpha*Omega[Lambdahat,:],np.concatenate((y, np.zeros(Lambdahat.size, 1))))));
+            temp = np.linalg.norm(y - np.dot(M,xhat), 2);
+            #disp(['fidelity error=', num2str(temp), ' lagmult=', num2str(lagmult)]);
+            if temp <= params['noise_level']:
+                was_feasible = 1;
+                if was_infeasible == 1:
+                    break;
+                else:
+                    lagmult = lagmult*lagmultfactor;
+            elif temp > params['noise_level']:
+                was_infeasible = 1;
+                if was_feasible == 1:
+                    xhat = xprev;
+                    break;
+                lagmult = lagmult/lagmultfactor;
+            if lagmult < lagmultmin or lagmult > lagmultmax:
+                break;
+            xprev = xhat;
+        arepr = np.dot(Omega[Lambdahat, :], xhat);
+        return xhat,arepr,lagmult;
+
+
+    ########################################################################
+    ## Computation using conjugate gradient method.
+    ########################################################################
+    #if strcmp(class(MH),'function_handle')
+    if hasattr(MH, '__call__'):
+        b = MH(y);
+    else:
+        b = np.dot(MH, y);
+
+    norm_b = np.linalg.norm(b, 2);
+    xhat = xinit;
+    xprev = xinit;
+    residual = TheHermitianMatrix(xhat, M, MH, Omega, OmegaH, Lambdahat, lagmult) - b;
+    direction = -residual;
+    iter = 0;
+
+    while iter < params.max_inner_iteration:
+        iter = iter + 1;
+        alpha = np.linalg.norm(residual,2)**2 / np.dot(direction.T, TheHermitianMatrix(direction, M, MH, Omega, OmegaH, Lambdahat, lagmult));
+        xhat = xhat + alpha*direction;
+        prev_residual = residual;
+        residual = TheHermitianMatrix(xhat, M, MH, Omega, OmegaH, Lambdahat, lagmult) - b;
+        beta = np.linalg.norm(residual,2)**2 / np.linalg.norm(prev_residual,2)**2;
+        direction = -residual + beta*direction;
+
+        if np.linalg.norm(residual,2)/norm_b < params['l2_accuracy']*(lagmult**(accuracy_adjustment_exponent)) or iter == params['max_inner_iteration']:
+            #if strcmp(class(M), 'function_handle')
+            if hasattr(M, '__call__'):
+                temp = np.linalg.norm(y-M(xhat), 2);
+            else:
+                temp = np.linalg.norm(y-np.dot(M,xhat), 2);
+
+            #if strcmp(class(Omega), 'function_handle')
+            if hasattr(Omega, '__call__'):
+                u = Omega(xhat);
+                u = math.sqrt(lagmult)*np.linalg.norm(u(Lambdahat), 2);
+            else:
+                u = math.sqrt(lagmult)*np.linalg.norm(Omega[Lambdahat,:]*xhat, 2);
+
+
+            #disp(['residual=', num2str(norm(residual,2)), ' norm_b=', num2str(norm_b), ' omegapart=', num2str(u), ' fidelity error=', num2str(temp), ' lagmult=', num2str(lagmult), ' iter=', num2str(iter)]);
+
+            if temp <= params['noise_level']:
+                was_feasible = 1;
+                if was_infeasible == 1:
+                    break;
+                else:
+                    lagmult = lagmultfactor*lagmult;
+                    residual = TheHermitianMatrix(xhat, M, MH, Omega, OmegaH, Lambdahat, lagmult) - b;
+                    direction = -residual;
+                    iter = 0;
+            elif temp > params['noise_level']:
+                lagmult = lagmult/lagmultfactor;
+                if was_feasible == 1:
+                    xhat = xprev;
+                    break;
+                was_infeasible = 1;
+                residual = TheHermitianMatrix(xhat, M, MH, Omega, OmegaH, Lambdahat, lagmult) - b;
+                direction = -residual;
+                iter = 0;
+            if lagmult > lagmultmax or lagmult < lagmultmin:
+                break;
+            xprev = xhat;
+        #elseif norm(xprev-xhat)/norm(xhat) < 1e-2
+        #    disp(['rel_change=', num2str(norm(xprev-xhat)/norm(xhat))]);
+        #    if strcmp(class(M), 'function_handle')
+        #        temp = norm(y-M(xhat), 2);
+        #    else
+        #        temp = norm(y-M*xhat, 2);
+        #    end
+    #
+    #        if temp > 1.2*params.noise_level
+    #            was_infeasible = 1;
+    #            lagmult = lagmult/lagmultfactor;
+    #            xprev = xhat;
+    #        end
+
+    #disp(['fidelity_error=', num2str(temp)]);
+    print 'fidelity_error=',temp
+    #if iter == params['max_inner_iteration']:
+        #disp('max_inner_iteration reached. l2_accuracy not achieved.');
+
+    ##
+    # Compute analysis representation for xhat
+    ##
+    #if strcmp(class(Omega),'function_handle')
+    if hasattr(Omega, '__call__'):
+        temp = Omega(xhat);
+        arepr = temp(Lambdahat);
+    else:    ## here Omega is assumed to be a matrix
+        arepr = np.dot(Omega[Lambdahat, :], xhat);
+
+    return xhat,arepr,lagmult
+
+
+##
+# This function computes (M'*M + lm*Omega(L,:)'*Omega(L,:)) * x.
+##
+#function w = TheHermitianMatrix(x, M, MH, Omega, OmegaH, L, lm)
+def TheHermitianMatrix(x, M, MH, Omega, OmegaH, L, lm):
+    #if strcmp(class(M), 'function_handle')
+    if hasattr(M, '__call__'):
+        w = MH(M(x));
+    else:    ## M and MH are matrices
+        w = np.dot(np.dot(MH, M), x);
+
+    if hasattr(Omega, '__call__'):
+        v = Omega(x);
+        vt = np.zeros(v.size);
+        vt[L] = v[L].copy();
+        w = w + lm*OmegaH(vt);
+    else:    ## Omega is assumed to be a matrix and OmegaH is its conjugate transpose
+        w = w + lm*np.dot(np.dot(OmegaH[:, L],Omega[L, :]),x);
+
+    return w
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/GAP_test.py	Sun Nov 06 20:58:11 2011 +0000
@@ -0,0 +1,69 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sun Nov 06 20:53:14 2011
+
+@author: Nic
+"""
+import numpy as np
+import numpy.linalg
+import scipy.io
+import unittest
+from pyCSalgos.SL0.SL0_approx import SL0_approx
+
+class SL0results(unittest.TestCase):
+  def testResults(self):
+    mdict = scipy.io.loadmat('SL0approxtestdata.mat')
+
+    # A = system matrix
+    # Y = matrix with measurements (on columns)
+    # sigmamin = vector with sigma_min
+    for A,Y,eps,sigmamin,Xr in zip(mdict['cellA'].squeeze(),mdict['cellY'].squeeze(),mdict['cellEps'].squeeze(),mdict['sigmamin'].squeeze(),mdict['cellXr'].squeeze()):
+      for i in np.arange(Y.shape[1]):
+
+        # Fix numpy error "LapackError: Parameter a has non-native byte order in lapack_lite.dgesdd"
+        A = A.newbyteorder('=')
+        Y = Y.newbyteorder('=')
+        eps = eps.newbyteorder('=')
+        sigmamin = sigmamin.newbyteorder('=')
+        Xr = Xr.newbyteorder('=')
+
+        xr = SL0_approx(A, Y[:,i], eps.squeeze()[i], sigmamin)
+
+        # check if found solution is the same as the correct cslution
+        diff = numpy.linalg.norm(xr - Xr[:,i])
+        self.assertTrue(diff < 1e-12)
+    #        err1 = numpy.linalg.norm(Y[:,i] - np.dot(A,xr))
+    #        err2 = numpy.linalg.norm(Y[:,i] - np.dot(A,Xr[:,i]))
+    #        norm1 = numpy.linalg.norm(xr,1)
+    #        norm2 = numpy.linalg.norm(Xr[:,i],1)
+    #
+    #        # Make a more robust condition:
+    #        #  OK;    if   solutions are close enough (diff < 1e-6)
+    #        #              or
+    #        #              (
+    #        #               Python solution fulfills the constraint better (or up to 1e-6 worse)
+    #        #                 and
+    #        #               Python solution has l1 norm no more than 1e-6 larger as the reference solution
+    #        #                 (i.e. either norm1 < norm2   or   norm1>norm2 not by more than 1e-6)
+    #        #              )
+    #        #
+    #        #  ERROR: else
+    #        differr  = err1 - err2    # intentionately no abs(), since err1` < err2 is good
+    #        diffnorm = norm1 - norm2  # intentionately no abs(), since norm1 < norm2 is good
+    #        if diff < 1e-6 or (differr < 1e-6 and (diffnorm < 1e-6)):
+    #          isok = True
+    #        else:
+    #          isok = False
+    #        self.assertTrue(isok)
+
+        #diff = numpy.linalg.norm(xr - Xr[:,i])
+        #if diff > 1e-6:
+        #    self.assertTrue(diff < 1e-6)
+
+
+if __name__ == "__main__":
+    #import cProfile
+    #cProfile.run('unittest.main()', 'profres')
+    unittest.main()
+    #suite = unittest.TestLoader().loadTestsFromTestCase(CompareResults)
+    #unittest.TextTestRunner(verbosity=2).run(suite)