Mercurial > hg > pycsalgos
changeset 67:a8d96e67717e
Added the Analysis-By-Synthesis algorithms used in the papers "Analysis-based sparse reconstruction with synthesis-based solvers", "Choosing Analysis or Synthesis Recovery for Sparse Reconstruction" and "A generalization of synthesis and analysis sparsity"
| author | Nic Cleju <nikcleju@gmail.com> |
|---|---|
| date | Tue, 09 Jul 2013 14:21:10 +0300 |
| parents | ee10ffb60428 |
| children | cab8a215f9a1 |
| files | pyCSalgos/ABS/ABSexact.py pyCSalgos/ABS/ABSlambda.py pyCSalgos/ABS/ABSmixed.py pyCSalgos/ABS/__init__.py pyCSalgos/ABS/algos.py pyCSalgos/GAP/GAP.py pyCSalgos/SL0/EllipseProj.py |
| diffstat | 6 files changed, 441 insertions(+), 3 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pyCSalgos/ABS/ABSexact.py Tue Jul 09 14:21:10 2013 +0300 @@ -0,0 +1,122 @@ +# -*- coding: utf-8 -*- +""" +Algorithms for exact analysis recovery based on synthesis solvers (a.k.a. Analysis by Synthesis, ABS). +Exact reconstruction. + +Author: Nicolae Cleju +""" +__author__ = "Nicolae Cleju" +__license__ = "GPL" +__email__ = "nikcleju@gmail.com" + + +import numpy + +# Import synthesis solvers from pyCSalgos package +import pyCSalgos.BP.l1eq_pd +import pyCSalgos.BP.cvxopt_lp +import pyCSalgos.OMP.omp_QR +import pyCSalgos.SL0.SL0 +import pyCSalgos.TST.RecommendedTST + + +def bp(y,M,Omega,x0, pdtol=1e-3, pdmaxiter=50, cgtol=1e-8, cgmaxiter=200, verbose=False): + """ + ABS-exact: Basis Pursuit (based on l1magic toolbox) + """ + N,n = Omega.shape + D = numpy.linalg.pinv(Omega) + U,S,Vt = numpy.linalg.svd(D) + Aextra = Vt[-(N-n):,:] + + # Create aggregate problem + Atilde = numpy.vstack((numpy.dot(M,D), Aextra)) + ytilde = numpy.concatenate((y,numpy.zeros(N-n))) + + return numpy.dot(D , pyCSalgos.BP.l1eq_pd.l1eq_pd(x0,Atilde,Atilde.T,ytilde, pdtol, pdmaxiter, cgtol, cgmaxiter, verbose)) + +def bp_cvxopt(y,M,Omega): + """ + ABS-exact: Basis Pursuit (based cvxopt) + """ + N,n = Omega.shape + D = numpy.linalg.pinv(Omega) + U,S,Vt = numpy.linalg.svd(D) + Aextra = Vt[-(N-n):,:] + + # Create aggregate problem + Atilde = numpy.vstack((numpy.dot(M,D), Aextra)) + ytilde = numpy.concatenate((y,numpy.zeros(N-n))) + + return numpy.dot(D , pyCSalgos.BP.cvxopt_lp.cvxopt_lp(ytilde, Atilde)) + + +def ompeps(y,M,Omega,epsilon): + """ + ABS-exact: OMP with stopping criterion residual < epsilon + """ + + N,n = Omega.shape + D = numpy.linalg.pinv(Omega) + U,S,Vt = numpy.linalg.svd(D) + Aextra = Vt[-(N-n):,:] + + # Create aggregate problem + Atilde = numpy.vstack((numpy.dot(M,D), Aextra)) + ytilde = numpy.concatenate((y,numpy.zeros(N-n))) + + opts = dict() + opts['stopCrit'] = 'mse' + opts['stopTol'] = epsilon + return numpy.dot(D , pyCSalgos.OMP.omp_QR.greed_omp_qr(ytilde,Atilde,Atilde.shape[1],opts)[0]) + +def ompk(y,M,Omega,k): + """ + ABS-exact: OMP with stopping criterion fixed number of atoms = k + """ + + N,n = Omega.shape + D = numpy.linalg.pinv(Omega) + U,S,Vt = numpy.linalg.svd(D) + Aextra = Vt[-(N-n):,:] + + # Create aggregate problem + Atilde = numpy.vstack((numpy.dot(M,D), Aextra)) + ytilde = numpy.concatenate((y,numpy.zeros(N-n))) + + opts = dict() + opts['stopTol'] = k + return numpy.dot(D , pyCSalgos.OMP.omp_QR.greed_omp_qr(ytilde,Atilde,Atilde.shape[1],opts)[0]) + +def sl0(y,M,Omega, sigma_min, sigma_decrease_factor=0.5, mu_0=2, L=3, true_s=None): + """ + ABS-exact: Smooth L0 (SL0) + """ + + N,n = Omega.shape + D = numpy.linalg.pinv(Omega) + U,S,Vt = numpy.linalg.svd(D) + Aextra = Vt[-(N-n):,:] + + # Create aggregate problem + Atilde = numpy.vstack((numpy.dot(M,D), Aextra)) + ytilde = numpy.concatenate((y,numpy.zeros(N-n))) + + return numpy.dot(D, pyCSalgos.SL0.SL0.SL0(Atilde,ytilde,sigma_min,sigma_decrease_factor,mu_0,L,true_s)) + +def tst_recom(y,M,Omega, nsweep=300, tol=0.00001, xinitial=None, ro=None): + """ + ABS-exact: Two Stage Thresholding (TST) with optimized parameters (see Maleki & Donoho) + """ + + N,n = Omega.shape + D = numpy.linalg.pinv(Omega) + U,S,Vt = numpy.linalg.svd(D) + Aextra = Vt[-(N-n):,:] + + # Create aggregate problem + Atilde = numpy.vstack((numpy.dot(M,D), Aextra)) + ytilde = numpy.concatenate((y,numpy.zeros(N-n))) + + return numpy.dot(D, pyCSalgos.TST.RecommendedTST.RecommendedTST(Atilde, ytilde, nsweep, tol, xinitial, ro)) + \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pyCSalgos/ABS/ABSlambda.py Tue Jul 09 14:21:10 2013 +0300 @@ -0,0 +1,81 @@ +# -*- coding: utf-8 -*- +""" +Algorithms for approximate analysis recovery based on synthesis solvers (a.k.a. Analysis by Synthesis, ABS). +Approximate reconstruction, ABS-lambda. + +Author: Nicolae Cleju +""" +__author__ = "Nicolae Cleju" +__license__ = "GPL" +__email__ = "nikcleju@gmail.com" + + +import numpy + +# Import synthesis solvers from pyCSalgos package +import pyCSalgos.BP.l1qc +import pyCSalgos.SL0.SL0_approx +import pyCSalgos.OMP.omp_QR +import pyCSalgos.TST.RecommendedTST + +def sl0(y,M,Omega,epsilon,lbd,sigma_min, sigma_decrease_factor=0.5, mu_0=2, L=3, A_pinv=None, true_s=None): + """ + ABS-lambda: Smooth L0 (SL0) + """ + N,n = Omega.shape + D = numpy.linalg.pinv(Omega) + U,S,Vt = numpy.linalg.svd(D) + aggDupper = numpy.dot(M,D) + aggDlower = Vt[-(N-n):,:] + aggD = numpy.vstack((aggDupper, lbd * aggDlower)) + aggy = numpy.concatenate((y, numpy.zeros(N-n))) + + #return numpy.dot(D, pyCSalgos.SL0.SL0_approx.SL0_approx(aggD,aggy,epsilon,sigma_min,sigma_decrease_factor,mu_0,L,A_pinv,true_s)) + return numpy.dot(D, pyCSalgos.SL0.SL0_approx.SL0_approx_dai(aggD,aggy,epsilon,sigma_min,sigma_decrease_factor,mu_0,L,A_pinv,true_s)) + +def bp(y,M,Omega,epsilon,lbd, x0, lbtol=1e-3, mu=10, cgtol=1e-8, cgmaxiter=200, verbose=False): + """ + ABS-lambda: Basis Pursuit (based on l1magic toolbox) + """ + N,n = Omega.shape + D = numpy.linalg.pinv(Omega) + U,S,Vt = numpy.linalg.svd(D) + aggDupper = numpy.dot(M,D) + aggDlower = Vt[-(N-n):,:] + aggD = numpy.vstack((aggDupper, lbd * aggDlower)) + aggy = numpy.concatenate((y, numpy.zeros(N-n))) + + return numpy.dot(D, pyCSalgos.BP.l1qc.l1qc_logbarrier(x0,aggD,aggD.T,aggy,epsilon, lbtol, mu, cgtol, cgmaxiter, verbose)) + +def ompeps(y,M,Omega,epsilon,lbd): + """ + ABS-lambda: OMP with stopping criterion residual < epsilon + """ + N,n = Omega.shape + D = numpy.linalg.pinv(Omega) + U,S,Vt = numpy.linalg.svd(D) + aggDupper = numpy.dot(M,D) + aggDlower = Vt[-(N-n):,:] + aggD = numpy.vstack((aggDupper, lbd * aggDlower)) + aggy = numpy.concatenate((y, numpy.zeros(N-n))) + + opts = dict() + opts['stopCrit'] = 'mse' + opts['stopTol'] = epsilon**2 / aggy.size + return numpy.dot(D, pyCSalgos.OMP.omp_QR.greed_omp_qr(aggy,aggD,aggD.shape[1],opts)[0]) + +def tst_recom(y,M,Omega,epsilon,lbd, nsweep=300, xinitial=None, ro=None): + """ + ABS-lambda: Two Stage Thresholding (TST) with optimized parameters (see Maleki & Donoho) + """ + N,n = Omega.shape + D = numpy.linalg.pinv(Omega) + U,S,Vt = numpy.linalg.svd(D) + aggDupper = numpy.dot(M,D) + aggDlower = Vt[-(N-n):,:] + aggD = numpy.vstack((aggDupper, lbd * aggDlower)) + aggy = numpy.concatenate((y, numpy.zeros(N-n))) + + tol = epsilon / numpy.linalg.norm(aggy) + return numpy.dot(D, pyCSalgos.TST.RecommendedTST.RecommendedTST(aggD, aggy, nsweep, tol, xinitial, ro)) + \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pyCSalgos/ABS/ABSmixed.py Tue Jul 09 14:21:10 2013 +0300 @@ -0,0 +1,45 @@ +# -*- coding: utf-8 -*- +""" +Algorithms for approximate analysis recovery based on synthesis solvers (a.k.a. Analysis by Synthesis, ABS). +Approximate reconstruction, ABS-mixed. + +Author: Nicolae Cleju +""" +__author__ = "Nicolae Cleju" +__license__ = "GPL" +__email__ = "nikcleju@gmail.com" + + +import numpy + +# Import synthesis solvers from pyCSalgos package +import pyCSalgos.BP.l1qec +import pyCSalgos.SL0.SL0_approx + +def bp(y,M,Omega,epsilon, x0, lbtol=1e-3, mu=10, cgtol=1e-8, cgmaxiter=200, verbose=False): + """ + ABS-mixed: Basis Pursuit (based on l1magic toolbox) + """ + N,n = Omega.shape + D = numpy.linalg.pinv(Omega) + U,S,Vt = numpy.linalg.svd(D) + Aeps = numpy.dot(M,D) + Aexact = Vt[-(N-n):,:] + + return numpy.dot(D , pyCSalgos.BP.l1qec.l1qec_logbarrier(x0,Aeps,Aeps.T,y,epsilon,Aexact,Aexact.T,numpy.zeros(N-n), lbtol, mu, cgtol, cgmaxiter, verbose)) + +def sl0(y,M,Omega,epsilon, sigma_min, sigma_decrease_factor=0.5, mu_0=2, L=3, Aeps_pinv=None, Aexact_pinv=None, true_s=None): + """ + ABS-mixed: Smooth L0 (SL0) + """ + N,n = Omega.shape + D = numpy.linalg.pinv(Omega) + U,S,Vt = numpy.linalg.svd(D) + Aeps = numpy.dot(M,D) + Aexact = Vt[-(N-n):,:] + + #return numpy.dot(D, pyCSalgos.SL0.SL0_approx.SL0_approx_analysis(Aeps,Aexact,y,epsilon,sigma_min,sigma_decrease_factor,mu_0,L,Aeps_pinv,Aexact_pinv,true_s)) + #return numpy.dot(D, pyCSalgos.SL0.SL0_approx.SL0_robust_analysis(Aeps,Aexact,y,epsilon,sigma_min,sigma_decrease_factor,mu_0,L,Aeps_pinv,Aexact_pinv,true_s)) + #return numpy.dot(D, pyCSalgos.SL0.SL0_approx.SL0_approx_analysis_unconstrained(Aeps,Aexact,y,epsilon,sigma_min,sigma_decrease_factor,mu_0,L,Aeps_pinv,Aexact_pinv,true_s)) + return numpy.dot(D, pyCSalgos.SL0.SL0_approx.SL0_approx_analysis_dai(Aeps,Aexact,y,epsilon,sigma_min,sigma_decrease_factor,mu_0,L,Aeps_pinv,Aexact_pinv,true_s)) + \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pyCSalgos/ABS/algos.py Tue Jul 09 14:21:10 2013 +0300 @@ -0,0 +1,188 @@ +# -*- coding: utf-8 -*- +""" +Define simple wrappers for algorithms, with similar header. +Specific algorithm parameters are defined inside here. + +Author: Nicolae Cleju +""" +__author__ = "Nicolae Cleju" +__license__ = "GPL" +__email__ = "nikcleju@gmail.com" + + +import numpy + +# Module with algorithms implemented in Python +import pyCSalgos +import pyCSalgos.GAP.GAP + +# Analysis by Synthesis - exact algorithms +import ABSexact +# Analysis by Synthesis - mixed algorithms +import ABSmixed +# Analysis by Synthesis - lambda algorithms +import ABSlambda + + +###--------------------------------- +### Exact reconstruction algorithms +###--------------------------------- +def run_exact_gap(y,M,Omega): + """ + Wrapper for GAP algorithm for exact analysis recovery + """ + gapparams = {"num_iteration" : 1000,\ + "greedy_level" : 0.9,\ + "stopping_coefficient_size" : 1e-4,\ + "l2solver" : 'pseudoinverse',\ + "noise_level": 1e-10} + return pyCSalgos.GAP.GAP.GAP(y,M,M.T,Omega,Omega.T,gapparams,numpy.zeros(Omega.shape[1]))[0] + +def run_exact_bp(y,M,Omega): + """ + Wrapper for BP algorithm for exact analysis recovery + Algorithm implementation is l1eq_pd() from l1-magic toolbox + """ + return ABSexact.bp(y,M,Omega,numpy.zeros(Omega.shape[0]), pdtol=1e-5, pdmaxiter = 100) + +def run_exact_bp_cvxopt(y,M,Omega): + """ + Wrapper for BP algorithm for exact analysis recovery + Algorithm implementation is using cvxopt linear programming + """ + return ABSexact.bp_cvxopt(y,M,Omega) + + +def run_exact_ompeps(y,M,Omega): + """ + Wrapper for OMP algorithm for exact analysis recovery, with stopping criterion = epsilon + """ + return ABSexact.ompeps(y,M,Omega,1e-9) + +#def run_exact_ompk(y,M,Omega) +# """ +# Wrapper for OMP algorithm for exact analysis recovery, with stopping criterion = fixed no. of atoms +# """ + +def run_exact_sl0(y,M,Omega): + """ + Wrapper for SL0 algorithm for exact analysis recovery + """ + sigma_min = 1e-12 + sigma_decrease_factor = 0.5 + mu_0 = 2 + L = 20 + return ABSexact.sl0(y,M,Omega, sigma_min, sigma_decrease_factor, mu_0, L) + +def run_exact_tst(y,M,Omega): + """ + Wrapper for TST algorithm (with default optimized params) for exact analysis recovery + """ + nsweep = 300 + tol = 1e-5 + return ABSexact.tst_recom(y,M,Omega, nsweep, tol) + + +###--------------------------------------- +### Approximate reconstruction algorithms +###--------------------------------------- +# 1. Native + +def run_gap(y,M,Omega,epsilon): + """ + Wrapper for GAP algorithm for approximate analysis recovery + """ + gapparams = {"num_iteration" : 1000,\ + "greedy_level" : 0.9,\ + "stopping_coefficient_size" : 1e-4,\ + "l2solver" : 'pseudoinverse',\ + "noise_level": epsilon} + return pyCSalgos.GAP.GAP.GAP(y,M,M.T,Omega,Omega.T,gapparams,numpy.zeros(Omega.shape[1]))[0] + +def run_nesta(y,M,Omega,epsilon): + """ + Wrapper for NESTA algorithm for approximate analysis recovery + """ + U,S,V = numpy.linalg.svd(M, full_matrices = True) + V = V.T # Make like Matlab + m,n = M.shape # Make like Matlab + S = numpy.hstack((numpy.diag(S), numpy.zeros((m,n-m)))) + + opt_muf = 1e-3 + optsUSV = {'U':U, 'S':S, 'V':V} + opts = {'U':Omega, 'Ut':Omega.T.copy(), 'USV':optsUSV, 'TolVar':1e-5, 'Verbose':0} + return pyCSalgos.NESTA.NESTA.NESTA(M, None, y, opt_muf, epsilon, opts)[0] + +# 2. ABS-mixed + +def run_mixed_sl0(y,M,Omega,epsilon): + """ + Wrapper for SL0-mixed algorithm for approximate analysis recovery + """ + sigma_min = 0.001 + sigma_decrease_factor = 0.5 + mu_0 = 2 + L = 10 + return ABSmixed.sl0(y,M,Omega,epsilon,sigma_min, sigma_decrease_factor, mu_0, L) + +def run_mixed_bp(y,M,Omega,epsilon): + """ + Wrapper for BP-mixed algorithm for approximate analysis recovery + """ + return ABSmixed.bp(y,M,Omega,epsilon, numpy.zeros(Omega.shape[0])) + +# 3. ABS-lambda + +def run_lambda_sl0(y,M,Omega,epsilon,lbd): + """ + Wrapper for SL0 algorithm within ABS-lambda approach for approximate analysis recovery + """ + sigma_min = 0.001 + sigma_decrease_factor = 0.5 + mu_0 = 2 + L = 10 + return ABSlambda.sl0(y,M,Omega,epsilon, lbd, sigma_min, sigma_decrease_factor, mu_0, L) + +def run_lambda_bp(y,M,Omega,epsilon,lbd): + """ + Wrapper for BP algorithm within ABS-lambda approach for approximate analysis recovery + """ + return ABSlambda.bp(y,M,Omega,epsilon,lbd,numpy.zeros(Omega.shape[0])) + +def run_lambda_ompeps(y,M,Omega,epsilon,lbd): + """ + Wrapper for OMP algorithm, with stopping criterion = epsilon, + for approximate analysis recovery within ABS-lambda approach + """ + return ABSlambda.ompeps(y,M,Omega,epsilon,lbd) + +def run_lambda_tst(y,M,Omega,epsilon,lbd): + """ + Wrapper for TST algorithm (with default optimized params) + for approximate analysis recovery within ABS-lambda approach + """ + nsweep = 300 + return ABSlambda.tst_recom(y,M,Omega,epsilon,lbd, nsweep) + + +### Define algorithm tuples: (function, name) +### Will be used in stdparams and in test scripts +## Exact recovery +exact_gap = (run_exact_gap, 'GAP') +exact_bp = (run_exact_bp, 'ABSexact_BP') +exact_bp_cvxopt = (run_exact_bp_cvxopt, 'ABSexact_BP_cvxopt') +exact_ompeps = (run_exact_ompeps, 'ABSexact_OMPeps') +exact_sl0 = (run_exact_sl0, 'ABSexact_SL0') +exact_tst = (run_exact_tst, 'ABSexact_TST') +## Approximate recovery +# Native +gap = (run_gap, 'GAP') +nesta = (run_nesta, 'NESTA') +# ABS-mixed +mixed_sl0 = (run_mixed_sl0, 'ABSmixed_SL0') +mixed_bp = (run_mixed_bp, 'ABSmixed_BP') +# ABS-lambda +lambda_sl0 = (run_lambda_sl0, 'ABSlambda_SL0') +lambda_bp = (run_lambda_bp, 'ABSlambda_BP') +lambda_ompeps = (run_lambda_ompeps, 'ABSlambda_OMPeps') +lambda_tst = (run_lambda_tst, 'ABSlambda_TST')
--- a/pyCSalgos/GAP/GAP.py Mon Apr 23 10:55:49 2012 +0300 +++ b/pyCSalgos/GAP/GAP.py Tue Jul 09 14:21:10 2013 +0300 @@ -9,6 +9,7 @@ import numpy import numpy.linalg import scipy as sp +import scipy.stats import math @@ -332,7 +333,8 @@ maxcoef = abscoef.max() if greedy_level >= 1: #qq = quantile(abscoef, 1.0-greedy_level/n); - qq = sp.stats.mstats.mquantile(abscoef, 1.0-greedy_level/n, 0.5, 0.5) + #qq = sp.stats.mstats.mquantiles(abscoef, 1.0-greedy_level/n, 0.5, 0.5) + qq = sp.stats.mstats.mquantiles(abscoef, 1.0-greedy_level/n) else: qq = maxcoef*greedy_level @@ -358,7 +360,7 @@ return 1 #if isfield(params, 'stopping_cosparsity') && length(Lambdahat) < params.stopping_cosparsity - if ('stopping_cosparsity' in params) and Lambdahat.size() < params['stopping_cosparsity']: + if ('stopping_cosparsity' in params) and Lambdahat.size < params['stopping_cosparsity']: return 1 return 0 \ No newline at end of file