Mercurial > hg > pycsalgos
changeset 0:8346c92b698f
Initial import
| author | nikcleju |
|---|---|
| date | Thu, 20 Oct 2011 19:36:24 +0000 |
| parents | |
| children | 2a2abf5092f8 |
| files | BP/cgsolve.m BP/l1qc.py BP/l1qc_logbarrier.m BP/l1qc_newton.m OMP/greed_omp_qr.m OMP/omp.py OMP/omp_downloaded.py OMP/omp_phase3.py |
| diffstat | 8 files changed, 3625 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/BP/cgsolve.m Thu Oct 20 19:36:24 2011 +0000 @@ -0,0 +1,76 @@ +% cgsolve.m +% +% Solve a symmetric positive definite system Ax = b via conjugate gradients. +% +% Usage: [x, res, iter] = cgsolve(A, b, tol, maxiter, verbose) +% +% A - Either an NxN matrix, or a function handle. +% +% b - N vector +% +% tol - Desired precision. Algorithm terminates when +% norm(Ax-b)/norm(b) < tol . +% +% maxiter - Maximum number of iterations. +% +% verbose - If 0, do not print out progress messages. +% If and integer greater than 0, print out progress every 'verbose' iters. +% +% Written by: Justin Romberg, Caltech +% Email: jrom@acm.caltech.edu +% Created: October 2005 +% + +function [x, res, iter] = cgsolve(A, b, tol, maxiter, verbose) + +if (nargin < 5), verbose = 1; end + +implicit = isa(A,'function_handle'); + +x = zeros(length(b),1); +r = b; +d = r; +delta = r'*r; +delta0 = b'*b; +numiter = 0; +bestx = x; +bestres = sqrt(delta/delta0); +while ((numiter < maxiter) && (delta > tol^2*delta0)) + + % q = A*d + if (implicit), q = A(d); else q = A*d; end + + alpha = delta/(d'*q); + x = x + alpha*d; + + if (mod(numiter+1,50) == 0) + % r = b - Aux*x + if (implicit), r = b - A(x); else r = b - A*x; end + else + r = r - alpha*q; + end + + deltaold = delta; + delta = r'*r; + beta = delta/deltaold; + d = r + beta*d; + numiter = numiter + 1; + if (sqrt(delta/delta0) < bestres) + bestx = x; + bestres = sqrt(delta/delta0); + end + + if ((verbose) && (mod(numiter,verbose)==0)) + disp(sprintf('cg: Iter = %d, Best residual = %8.3e, Current residual = %8.3e', ... + numiter, bestres, sqrt(delta/delta0))); + end + +end + +if (verbose) + disp(sprintf('cg: Iterations = %d, best residual = %14.8e', numiter, bestres)); +end +x = bestx; +res = bestres; +iter = numiter; +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/BP/l1qc.py Thu Oct 20 19:36:24 2011 +0000 @@ -0,0 +1,705 @@ + +import numpy as np +import scipy.linalg +import math + + +#function [x, res, iter] = cgsolve(A, b, tol, maxiter, verbose) +def cgsolve(A, b, tol, maxiter, verbose=1): + # Solve a symmetric positive definite system Ax = b via conjugate gradients. + # + # Usage: [x, res, iter] = cgsolve(A, b, tol, maxiter, verbose) + # + # A - Either an NxN matrix, or a function handle. + # + # b - N vector + # + # tol - Desired precision. Algorithm terminates when + # norm(Ax-b)/norm(b) < tol . + # + # maxiter - Maximum number of iterations. + # + # verbose - If 0, do not print out progress messages. + # If and integer greater than 0, print out progress every 'verbose' iters. + # + # Written by: Justin Romberg, Caltech + # Email: jrom@acm.caltech.edu + # Created: October 2005 + # + + #--------------------- + # Original Matab code: + # + #if (nargin < 5), verbose = 1; end + # + #implicit = isa(A,'function_handle'); + # + #x = zeros(length(b),1); + #r = b; + #d = r; + #delta = r'*r; + #delta0 = b'*b; + #numiter = 0; + #bestx = x; + #bestres = sqrt(delta/delta0); + #while ((numiter < maxiter) && (delta > tol^2*delta0)) + # + # # q = A*d + # if (implicit), q = A(d); else q = A*d; end + # + # alpha = delta/(d'*q); + # x = x + alpha*d; + # + # if (mod(numiter+1,50) == 0) + # # r = b - Aux*x + # if (implicit), r = b - A(x); else r = b - A*x; end + # else + # r = r - alpha*q; + # end + # + # deltaold = delta; + # delta = r'*r; + # beta = delta/deltaold; + # d = r + beta*d; + # numiter = numiter + 1; + # if (sqrt(delta/delta0) < bestres) + # bestx = x; + # bestres = sqrt(delta/delta0); + # end + # + # if ((verbose) && (mod(numiter,verbose)==0)) + # disp(sprintf('cg: Iter = #d, Best residual = #8.3e, Current residual = #8.3e', ... + # numiter, bestres, sqrt(delta/delta0))); + # end + # + #end + # + #if (verbose) + # disp(sprintf('cg: Iterations = #d, best residual = #14.8e', numiter, bestres)); + #end + #x = bestx; + #res = bestres; + #iter = numiter; + + # End of original Matab code + #---------------------------- + + #if (nargin < 5), verbose = 1; end + # Optional argument + + #implicit = isa(A,'function_handle'); + if hasattr(A, '__call__'): + implicit = True + else: + implicit = False + + x = np.zeros(b.size) + r = b.copy() + d = r.copy() + delta = np.vdot(r,r) + delta0 = np.vdot(b,b) + numiter = 0 + bestx = x.copy() + bestres = math.sqrt(delta/delta0) + while (numiter < maxiter) and (delta > tol**2*delta0): + + # q = A*d + #if (implicit), q = A(d); else q = A*d; end + if implicit: + q = A(d) + else: + q = np.dot(A,d) + + alpha = delta/np.vdot(d,q) + x = x + alpha*d + + if divmod(numiter+1,50)[1] == 0: + # r = b - Aux*x + #if (implicit), r = b - A(x); else r = b - A*x; end + if implicit: + r = b - A(x) + else: + r = b - np.dot(A,x) + else: + r = r - alpha*q + #end + + deltaold = delta; + delta = np.vdot(r,r) + beta = delta/deltaold; + d = r + beta*d + numiter = numiter + 1 + if (math.sqrt(delta/delta0) < bestres): + bestx = x.copy() + bestres = math.sqrt(delta/delta0) + #end + + if ((verbose) and (divmod(numiter,verbose)[1]==0)): + #disp(sprintf('cg: Iter = #d, Best residual = #8.3e, Current residual = #8.3e', ... + # numiter, bestres, sqrt(delta/delta0))); + print 'cg: Iter = ',numiter,', Best residual = ',bestres,', Current residual = ',math.sqrt(delta/delta0) + #end + + #end + + if (verbose): + #disp(sprintf('cg: Iterations = #d, best residual = #14.8e', numiter, bestres)); + print 'cg: Iterations = ',numiter,', best residual = ',bestres + #end + x = bestx.copy() + res = bestres + iter = numiter + + return x,res,iter + + + +#function [xp, up, niter] = l1qc_newton(x0, u0, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter) +def l1qc_newton(x0, u0, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter): + # Newton algorithm for log-barrier subproblems for l1 minimization + # with quadratic constraints. + # + # Usage: + # [xp,up,niter] = l1qc_newton(x0, u0, A, At, b, epsilon, tau, + # newtontol, newtonmaxiter, cgtol, cgmaxiter) + # + # x0,u0 - starting points + # + # A - Either a handle to a function that takes a N vector and returns a K + # vector , or a KxN matrix. If A is a function handle, the algorithm + # operates in "largescale" mode, solving the Newton systems via the + # Conjugate Gradients algorithm. + # + # At - Handle to a function that takes a K vector and returns an N vector. + # If A is a KxN matrix, At is ignored. + # + # b - Kx1 vector of observations. + # + # epsilon - scalar, constraint relaxation parameter + # + # tau - Log barrier parameter. + # + # newtontol - Terminate when the Newton decrement is <= newtontol. + # Default = 1e-3. + # + # newtonmaxiter - Maximum number of iterations. + # Default = 50. + # + # cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix. + # Default = 1e-8. + # + # cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored + # if A is a matrix. + # Default = 200. + # + # Written by: Justin Romberg, Caltech + # Email: jrom@acm.caltech.edu + # Created: October 2005 + # + + #--------------------- + # Original Matab code: + + ## check if the matrix A is implicit or explicit + #largescale = isa(A,'function_handle'); + # + ## line search parameters + #alpha = 0.01; + #beta = 0.5; + # + #if (~largescale), AtA = A'*A; end + # + ## initial point + #x = x0; + #u = u0; + #if (largescale), r = A(x) - b; else r = A*x - b; end + #fu1 = x - u; + #fu2 = -x - u; + #fe = 1/2*(r'*r - epsilon^2); + #f = sum(u) - (1/tau)*(sum(log(-fu1)) + sum(log(-fu2)) + log(-fe)); + # + #niter = 0; + #done = 0; + #while (~done) + # + # if (largescale), atr = At(r); else atr = A'*r; end + # + # ntgz = 1./fu1 - 1./fu2 + 1/fe*atr; + # ntgu = -tau - 1./fu1 - 1./fu2; + # gradf = -(1/tau)*[ntgz; ntgu]; + # + # sig11 = 1./fu1.^2 + 1./fu2.^2; + # sig12 = -1./fu1.^2 + 1./fu2.^2; + # sigx = sig11 - sig12.^2./sig11; + # + # w1p = ntgz - sig12./sig11.*ntgu; + # if (largescale) + # h11pfun = @(z) sigx.*z - (1/fe)*At(A(z)) + 1/fe^2*(atr'*z)*atr; + # [dx, cgres, cgiter] = cgsolve(h11pfun, w1p, cgtol, cgmaxiter, 0); + # if (cgres > 1/2) + # disp('Cannot solve system. Returning previous iterate. (See Section 4 of notes for more information.)'); + # xp = x; up = u; + # return + # end + # Adx = A(dx); + # else + # H11p = diag(sigx) - (1/fe)*AtA + (1/fe)^2*atr*atr'; + # opts.POSDEF = true; opts.SYM = true; + # [dx,hcond] = linsolve(H11p, w1p, opts); + # if (hcond < 1e-14) + # disp('Matrix ill-conditioned. Returning previous iterate. (See Section 4 of notes for more information.)'); + # xp = x; up = u; + # return + # end + # Adx = A*dx; + # end + # du = (1./sig11).*ntgu - (sig12./sig11).*dx; + # + # # minimum step size that stays in the interior + # ifu1 = find((dx-du) > 0); ifu2 = find((-dx-du) > 0); + # aqe = Adx'*Adx; bqe = 2*r'*Adx; cqe = r'*r - epsilon^2; + # smax = min(1,min([... + # -fu1(ifu1)./(dx(ifu1)-du(ifu1)); -fu2(ifu2)./(-dx(ifu2)-du(ifu2)); ... + # (-bqe+sqrt(bqe^2-4*aqe*cqe))/(2*aqe) + # ])); + # s = (0.99)*smax; + # + # # backtracking line search + # suffdec = 0; + # backiter = 0; + # while (~suffdec) + # xp = x + s*dx; up = u + s*du; rp = r + s*Adx; + # fu1p = xp - up; fu2p = -xp - up; fep = 1/2*(rp'*rp - epsilon^2); + # fp = sum(up) - (1/tau)*(sum(log(-fu1p)) + sum(log(-fu2p)) + log(-fep)); + # flin = f + alpha*s*(gradf'*[dx; du]); + # suffdec = (fp <= flin); + # s = beta*s; + # backiter = backiter + 1; + # if (backiter > 32) + # disp('Stuck on backtracking line search, returning previous iterate. (See Section 4 of notes for more information.)'); + # xp = x; up = u; + # return + # end + # end + # + # # set up for next iteration + # x = xp; u = up; r = rp; + # fu1 = fu1p; fu2 = fu2p; fe = fep; f = fp; + # + # lambda2 = -(gradf'*[dx; du]); + # stepsize = s*norm([dx; du]); + # niter = niter + 1; + # done = (lambda2/2 < newtontol) | (niter >= newtonmaxiter); + # + # disp(sprintf('Newton iter = #d, Functional = #8.3f, Newton decrement = #8.3f, Stepsize = #8.3e', ... + # niter, f, lambda2/2, stepsize)); + # if (largescale) + # disp(sprintf(' CG Res = #8.3e, CG Iter = #d', cgres, cgiter)); + # else + # disp(sprintf(' H11p condition number = #8.3e', hcond)); + # end + # + #end + + # End of original Matab code + #---------------------------- + + # check if the matrix A is implicit or explicit + #largescale = isa(A,'function_handle'); + if hasattr(A, '__call__'): + largescale = True + else: + largescale = False + + # line search parameters + alpha = 0.01 + beta = 0.5 + + #if (~largescale), AtA = A'*A; end + if not largescale: + AtA = np.dot(A.T,A) + + # initial point + x = x0.copy() + u = u0.copy() + #if (largescale), r = A(x) - b; else r = A*x - b; end + if largescale: + r = A(x) - b + else: + r = np.dot(A,x) - b + + fu1 = x - u + fu2 = -x - u + fe = 1/2*(np.vdot(r,r) - epsilon**2) + f = u.sum() - (1/tau)*(np.log(-fu1).sum() + np.log(-fu2).sum() + math.log(-fe)) + + niter = 0 + done = 0 + while not done: + + #if (largescale), atr = At(r); else atr = A'*r; end + if largescale: + atr = At(r) + else: + atr = np.dot(A.T,r) + + #ntgz = 1./fu1 - 1./fu2 + 1/fe*atr; + ntgz = 1.0/fu1 - 1.0/fu2 + 1.0/fe*atr + #ntgu = -tau - 1./fu1 - 1./fu2; + ntgu = -tau - 1.0/fu1 - 1.0/fu2 + #gradf = -(1/tau)*[ntgz; ntgu]; + gradf = -(1.0/tau)*np.concatenate((ntgz, ntgu),0) + + #sig11 = 1./fu1.^2 + 1./fu2.^2; + sig11 = 1.0/(fu1**2) + 1.0/(fu2**2) + #sig12 = -1./fu1.^2 + 1./fu2.^2; + sig12 = -1.0/(fu1**2) + 1.0/(fu2**2) + #sigx = sig11 - sig12.^2./sig11; + sigx = sig11 - (sig12**2)/sig11 + + #w1p = ntgz - sig12./sig11.*ntgu; + w1p = ntgz - sig12/sig11*ntgu + if largescale: + #h11pfun = @(z) sigx.*z - (1/fe)*At(A(z)) + 1/fe^2*(atr'*z)*atr; + h11pfun = lambda z: sigx*z - (1.0/fe)*At(A(z)) + 1.0/(fe**2)*np.dot(np.dot(atr.T,z),atr) + dx,cgres,cgiter = cgsolve(h11pfun, w1p, cgtol, cgmaxiter, 0) + if (cgres > 1/2): + print 'Cannot solve system. Returning previous iterate. (See Section 4 of notes for more information.)' + xp = x.copy() + up = u.copy() + return xp,up,niter + #end + Adx = A(dx) + else: + #H11p = diag(sigx) - (1/fe)*AtA + (1/fe)^2*atr*atr'; + H11p = np.diag(sigx) - (1/fe)*AtA + (1/fe)**2*np.dot(atr,atr.T) + #opts.POSDEF = true; opts.SYM = true; + #[dx,hcond] = linsolve(H11p, w1p, opts); + #if (hcond < 1e-14) + # disp('Matrix ill-conditioned. Returning previous iterate. (See Section 4 of notes for more information.)'); + # xp = x; up = u; + # return + #end + try: + dx = scipy.linalg.solve(H11p, w1p, sym_pos=True) + hcond = 1.0/scipy.linalg.cond(H11p) + except scipy.linalg.LinAlgError: + print 'Matrix ill-conditioned. Returning previous iterate. (See Section 4 of notes for more information.)' + xp = x.copy() + up = u.copy() + return xp,up,niter + if hcond < 1e-14: + print 'Matrix ill-conditioned. Returning previous iterate. (See Section 4 of notes for more information.)' + xp = x.copy() + up = u.copy() + return xp,up,niter + + #Adx = A*dx; + Adx = np.dot(A,dx) + #end + #du = (1./sig11).*ntgu - (sig12./sig11).*dx; + du = (1.0/sig11)*ntgu - (sig12/sig11)*dx; + + # minimum step size that stays in the interior + #ifu1 = find((dx-du) > 0); ifu2 = find((-dx-du) > 0); + ifu1 = np.nonzero((dx-du)>0) + ifu2 = np.nonzero((-dx-du)>0) + #aqe = Adx'*Adx; bqe = 2*r'*Adx; cqe = r'*r - epsilon^2; + aqe = np.dot(Adx.T,Adx) + bqe = 2*np.dot(r.T,Adx) + cqe = np.vdot(r,r) - epsilon**2 + #smax = min(1,min([... + # -fu1(ifu1)./(dx(ifu1)-du(ifu1)); -fu2(ifu2)./(-dx(ifu2)-du(ifu2)); ... + # (-bqe+sqrt(bqe^2-4*aqe*cqe))/(2*aqe) + # ])); + smax = min(1,np.concatenate((-fu1[ifu1]/(dx[ifu1]-du[ifu1]) , -fu2[ifu2]/(-dx[ifu2]-du[ifu2]) , (-bqe + math.sqrt(bqe**2-4*aqe*cqe))/(2*aqe)),0).min()) + + s = (0.99)*smax + + # backtracking line search + suffdec = 0 + backiter = 0 + while not suffdec: + #xp = x + s*dx; up = u + s*du; rp = r + s*Adx; + xp = x + s*dx + up = u + s*du + rp = r + s*Adx + #fu1p = xp - up; fu2p = -xp - up; fep = 1/2*(rp'*rp - epsilon^2); + fu1p = xp - up + fu2p = -xp - up + fep = 0.5*(np.vdot(r,r) - epsilon**2) + #fp = sum(up) - (1/tau)*(sum(log(-fu1p)) + sum(log(-fu2p)) + log(-fep)); + fp = up.sum() - (1.0/tau)*(np.log(-fu1p).sum() + np.log(-fu2p).sum() + log(-fep)) + #flin = f + alpha*s*(gradf'*[dx; du]); + flin = f + alpha*s*np.dot(gradf.T , np.concatenate((dx,du),0)) + #suffdec = (fp <= flin); + if fp <= flin: + suffdec = True + else: + suffdec = False + + s = beta*s + backiter = backiter + 1 + if (backiter > 32): + print 'Stuck on backtracking line search, returning previous iterate. (See Section 4 of notes for more information.)' + xp = x.copy() + up = u.copy() + return xp,up,niter + #end + #end + + # set up for next iteration + #x = xp; u = up; r = rp; + x = xp.copy() + u = up.copy() + r = rp.copy() + #fu1 = fu1p; fu2 = fu2p; fe = fep; f = fp; + fu1 = fu1p.copy() + fu2 = fu2p.copy() + fe = fep + f = fp + + #lambda2 = -(gradf'*[dx; du]); + lambda2 = -np.dot(gradf , np.concatenate((dx,du),0)) + #stepsize = s*norm([dx; du]); + stepsize = s * np.linalg.norm(np.concatenate((dx,du),0)) + niter = niter + 1 + #done = (lambda2/2 < newtontol) | (niter >= newtonmaxiter); + if lambda2/2.0 < newtontol or niter >= newtonmaxiter: + done = 1 + else: + done = 0 + + #disp(sprintf('Newton iter = #d, Functional = #8.3f, Newton decrement = #8.3f, Stepsize = #8.3e', ... + print 'Newton iter = ',niter,', Functional = ',f,', Newton decrement = ',lambda2/2.0,', Stepsize = ',stepsize + + if largescale: + print ' CG Res = ',cgres,', CG Iter = ',cgiter + else: + print ' H11p condition number = ',hcond + #end + + #end + return xp,up,niter + +#function xp = l1qc_logbarrier(x0, A, At, b, epsilon, lbtol, mu, cgtol, cgmaxiter) +def l1qc_logbarrier(x0, A, At, b, epsilon, lbtol=1e-3, mu=10, cgtol=1e-8, cgmaxiter=200): + # Solve quadratically constrained l1 minimization: + # min ||x||_1 s.t. ||Ax - b||_2 <= \epsilon + # + # Reformulate as the second-order cone program + # min_{x,u} sum(u) s.t. x - u <= 0, + # -x - u <= 0, + # 1/2(||Ax-b||^2 - \epsilon^2) <= 0 + # and use a log barrier algorithm. + # + # Usage: xp = l1qc_logbarrier(x0, A, At, b, epsilon, lbtol, mu, cgtol, cgmaxiter) + # + # x0 - Nx1 vector, initial point. + # + # A - Either a handle to a function that takes a N vector and returns a K + # vector , or a KxN matrix. If A is a function handle, the algorithm + # operates in "largescale" mode, solving the Newton systems via the + # Conjugate Gradients algorithm. + # + # At - Handle to a function that takes a K vector and returns an N vector. + # If A is a KxN matrix, At is ignored. + # + # b - Kx1 vector of observations. + # + # epsilon - scalar, constraint relaxation parameter + # + # lbtol - The log barrier algorithm terminates when the duality gap <= lbtol. + # Also, the number of log barrier iterations is completely + # determined by lbtol. + # Default = 1e-3. + # + # mu - Factor by which to increase the barrier constant at each iteration. + # Default = 10. + # + # cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix. + # Default = 1e-8. + # + # cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored + # if A is a matrix. + # Default = 200. + # + # Written by: Justin Romberg, Caltech + # Email: jrom@acm.caltech.edu + # Created: October 2005 + # + + #--------------------- + # Original Matab code: + + #largescale = isa(A,'function_handle'); + # + #if (nargin < 6), lbtol = 1e-3; end + #if (nargin < 7), mu = 10; end + #if (nargin < 8), cgtol = 1e-8; end + #if (nargin < 9), cgmaxiter = 200; end + # + #newtontol = lbtol; + #newtonmaxiter = 50; + # + #N = length(x0); + # + ## starting point --- make sure that it is feasible + #if (largescale) + # if (norm(A(x0)-b) > epsilon) + # disp('Starting point infeasible; using x0 = At*inv(AAt)*y.'); + # AAt = @(z) A(At(z)); + # w = cgsolve(AAt, b, cgtol, cgmaxiter, 0); + # if (cgres > 1/2) + # disp('A*At is ill-conditioned: cannot find starting point'); + # xp = x0; + # return; + # end + # x0 = At(w); + # end + #else + # if (norm(A*x0-b) > epsilon) + # disp('Starting point infeasible; using x0 = At*inv(AAt)*y.'); + # opts.POSDEF = true; opts.SYM = true; + # [w, hcond] = linsolve(A*A', b, opts); + # if (hcond < 1e-14) + # disp('A*At is ill-conditioned: cannot find starting point'); + # xp = x0; + # return; + # end + # x0 = A'*w; + # end + #end + #x = x0; + #u = (0.95)*abs(x0) + (0.10)*max(abs(x0)); + # + #disp(sprintf('Original l1 norm = #.3f, original functional = #.3f', sum(abs(x0)), sum(u))); + # + ## choose initial value of tau so that the duality gap after the first + ## step will be about the origial norm + #tau = max((2*N+1)/sum(abs(x0)), 1); + # + #lbiter = ceil((log(2*N+1)-log(lbtol)-log(tau))/log(mu)); + #disp(sprintf('Number of log barrier iterations = #d\n', lbiter)); + # + #totaliter = 0; + # + ## Added by Nic + #if lbiter == 0 + # xp = zeros(size(x0)); + #end + # + #for ii = 1:lbiter + # + # [xp, up, ntiter] = l1qc_newton(x, u, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter); + # totaliter = totaliter + ntiter; + # + # disp(sprintf('\nLog barrier iter = #d, l1 = #.3f, functional = #8.3f, tau = #8.3e, total newton iter = #d\n', ... + # ii, sum(abs(xp)), sum(up), tau, totaliter)); + # + # x = xp; + # u = up; + # + # tau = mu*tau; + # + #end + # + # End of original Matab code + #---------------------------- + + #largescale = isa(A,'function_handle'); + if hasattr(A, '__call__'): + largescale = True + else: + largescale = False + + # if (nargin < 6), lbtol = 1e-3; end + # if (nargin < 7), mu = 10; end + # if (nargin < 8), cgtol = 1e-8; end + # if (nargin < 9), cgmaxiter = 200; end + # Nic: added them as optional parameteres + + newtontol = lbtol + newtonmaxiter = 50 + + #N = length(x0); + N = x0.size() + + # starting point --- make sure that it is feasible + if largescale: + if np.linalg.norm(A(x0) - b) > epsilon: + print 'Starting point infeasible; using x0 = At*inv(AAt)*y.' + #AAt = @(z) A(At(z)); + AAt = lambda z: A(At(z)) + # TODO: implement cgsolve + w,cgres,cgiter = cgsolve(AAt, b, cgtol, cgmaxiter, 0) + if (cgres > 1/2): + print 'A*At is ill-conditioned: cannot find starting point' + xp = x0.copy() + return xp + #end + x0 = At(w) + #end + else: + if np.linalg.norm( np.dot(A,x0) - b ) > epsilon: + print 'Starting point infeasible; using x0 = At*inv(AAt)*y.' + #opts.POSDEF = true; opts.SYM = true; + #[w, hcond] = linsolve(A*A', b, opts); + #if (hcond < 1e-14) + # disp('A*At is ill-conditioned: cannot find starting point'); + # xp = x0; + # return; + #end + try: + w = scipy.linalg.solve(np.dot(A,A.T), b, sym_pos=True) + hcond = 1.0/scipy.linalg.cond(np.dot(A,A.T)) + except scipy.linalg.LinAlgError: + print 'A*At is ill-conditioned: cannot find starting point' + xp = x0.copy() + return xp + if hcond < 1e-14: + print 'A*At is ill-conditioned: cannot find starting point' + xp = x0.copy() + return xp + #x0 = A'*w; + x0 = np.dot(A.T, w) + #end + #end + x = x0.copy() + u = (0.95)*np.abs(x0) + (0.10)*np.abs(x0).max() + + #disp(sprintf('Original l1 norm = #.3f, original functional = #.3f', sum(abs(x0)), sum(u))); + print 'Original l1 norm = ',np.abs(x0).sum(),'original functional = ',u.sum() + + # choose initial value of tau so that the duality gap after the first + # step will be about the origial norm + tau = max(((2*N+1)/np.abs(x0).sum()), 1) + + lbiter = math.ceil((math.log(2*N+1)-math.log(lbtol)-math.log(tau))/math.log(mu)) + #disp(sprintf('Number of log barrier iterations = #d\n', lbiter)); + print 'Number of log barrier iterations = ',lbiter + + totaliter = 0 + + # Added by Nic, to fix some crashing + if lbiter == 0: + xp = np.zeros(x0.size) + #end + + #for ii = 1:lbiter + for ii in np.arange(lbiter): + + xp,up,ntiter = l1qc_newton(x, u, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter) + totaliter = totaliter + ntiter + + #disp(sprintf('\nLog barrier iter = #d, l1 = #.3f, functional = #8.3f, tau = #8.3e, total newton iter = #d\n', ... + # ii, sum(abs(xp)), sum(up), tau, totaliter)); + print 'Log barrier iter = ',ii,', l1 = ',np.abs(xp).sum(),', functional = ',up.sum(),', tau = ',tau,', total newton iter = ',totaliter + x = xp.copy() + u = up.copy() + + tau = mu*tau + + #end + return xp +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/BP/l1qc_logbarrier.m Thu Oct 20 19:36:24 2011 +0000 @@ -0,0 +1,121 @@ +% l1qc_logbarrier.m +% +% Solve quadratically constrained l1 minimization: +% min ||x||_1 s.t. ||Ax - b||_2 <= \epsilon +% +% Reformulate as the second-order cone program +% min_{x,u} sum(u) s.t. x - u <= 0, +% -x - u <= 0, +% 1/2(||Ax-b||^2 - \epsilon^2) <= 0 +% and use a log barrier algorithm. +% +% Usage: xp = l1qc_logbarrier(x0, A, At, b, epsilon, lbtol, mu, cgtol, cgmaxiter) +% +% x0 - Nx1 vector, initial point. +% +% A - Either a handle to a function that takes a N vector and returns a K +% vector , or a KxN matrix. If A is a function handle, the algorithm +% operates in "largescale" mode, solving the Newton systems via the +% Conjugate Gradients algorithm. +% +% At - Handle to a function that takes a K vector and returns an N vector. +% If A is a KxN matrix, At is ignored. +% +% b - Kx1 vector of observations. +% +% epsilon - scalar, constraint relaxation parameter +% +% lbtol - The log barrier algorithm terminates when the duality gap <= lbtol. +% Also, the number of log barrier iterations is completely +% determined by lbtol. +% Default = 1e-3. +% +% mu - Factor by which to increase the barrier constant at each iteration. +% Default = 10. +% +% cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix. +% Default = 1e-8. +% +% cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored +% if A is a matrix. +% Default = 200. +% +% Written by: Justin Romberg, Caltech +% Email: jrom@acm.caltech.edu +% Created: October 2005 +% + +function xp = l1qc_logbarrier(x0, A, At, b, epsilon, lbtol, mu, cgtol, cgmaxiter) + +largescale = isa(A,'function_handle'); + +if (nargin < 6), lbtol = 1e-3; end +if (nargin < 7), mu = 10; end +if (nargin < 8), cgtol = 1e-8; end +if (nargin < 9), cgmaxiter = 200; end + +newtontol = lbtol; +newtonmaxiter = 50; + +N = length(x0); + +% starting point --- make sure that it is feasible +if (largescale) + if (norm(A(x0)-b) > epsilon) + disp('Starting point infeasible; using x0 = At*inv(AAt)*y.'); + AAt = @(z) A(At(z)); + w = cgsolve(AAt, b, cgtol, cgmaxiter, 0); + if (cgres > 1/2) + disp('A*At is ill-conditioned: cannot find starting point'); + xp = x0; + return; + end + x0 = At(w); + end +else + if (norm(A*x0-b) > epsilon) + disp('Starting point infeasible; using x0 = At*inv(AAt)*y.'); + opts.POSDEF = true; opts.SYM = true; + [w, hcond] = linsolve(A*A', b, opts); + if (hcond < 1e-14) + disp('A*At is ill-conditioned: cannot find starting point'); + xp = x0; + return; + end + x0 = A'*w; + end +end +x = x0; +u = (0.95)*abs(x0) + (0.10)*max(abs(x0)); + +disp(sprintf('Original l1 norm = %.3f, original functional = %.3f', sum(abs(x0)), sum(u))); + +% choose initial value of tau so that the duality gap after the first +% step will be about the origial norm +tau = max((2*N+1)/sum(abs(x0)), 1); + +lbiter = ceil((log(2*N+1)-log(lbtol)-log(tau))/log(mu)); +disp(sprintf('Number of log barrier iterations = %d\n', lbiter)); + +totaliter = 0; + +% Added by Nic +if lbiter == 0 + xp = zeros(size(x0)); +end + +for ii = 1:lbiter + + [xp, up, ntiter] = l1qc_newton(x, u, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter); + totaliter = totaliter + ntiter; + + disp(sprintf('\nLog barrier iter = %d, l1 = %.3f, functional = %8.3f, tau = %8.3e, total newton iter = %d\n', ... + ii, sum(abs(xp)), sum(up), tau, totaliter)); + + x = xp; + u = up; + + tau = mu*tau; + +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/BP/l1qc_newton.m Thu Oct 20 19:36:24 2011 +0000 @@ -0,0 +1,150 @@ +% l1qc_newton.m +% +% Newton algorithm for log-barrier subproblems for l1 minimization +% with quadratic constraints. +% +% Usage: +% [xp,up,niter] = l1qc_newton(x0, u0, A, At, b, epsilon, tau, +% newtontol, newtonmaxiter, cgtol, cgmaxiter) +% +% x0,u0 - starting points +% +% A - Either a handle to a function that takes a N vector and returns a K +% vector , or a KxN matrix. If A is a function handle, the algorithm +% operates in "largescale" mode, solving the Newton systems via the +% Conjugate Gradients algorithm. +% +% At - Handle to a function that takes a K vector and returns an N vector. +% If A is a KxN matrix, At is ignored. +% +% b - Kx1 vector of observations. +% +% epsilon - scalar, constraint relaxation parameter +% +% tau - Log barrier parameter. +% +% newtontol - Terminate when the Newton decrement is <= newtontol. +% Default = 1e-3. +% +% newtonmaxiter - Maximum number of iterations. +% Default = 50. +% +% cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix. +% Default = 1e-8. +% +% cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored +% if A is a matrix. +% Default = 200. +% +% Written by: Justin Romberg, Caltech +% Email: jrom@acm.caltech.edu +% Created: October 2005 +% + + +function [xp, up, niter] = l1qc_newton(x0, u0, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter) + +% check if the matrix A is implicit or explicit +largescale = isa(A,'function_handle'); + +% line search parameters +alpha = 0.01; +beta = 0.5; + +if (~largescale), AtA = A'*A; end + +% initial point +x = x0; +u = u0; +if (largescale), r = A(x) - b; else r = A*x - b; end +fu1 = x - u; +fu2 = -x - u; +fe = 1/2*(r'*r - epsilon^2); +f = sum(u) - (1/tau)*(sum(log(-fu1)) + sum(log(-fu2)) + log(-fe)); + +niter = 0; +done = 0; +while (~done) + + if (largescale), atr = At(r); else atr = A'*r; end + + ntgz = 1./fu1 - 1./fu2 + 1/fe*atr; + ntgu = -tau - 1./fu1 - 1./fu2; + gradf = -(1/tau)*[ntgz; ntgu]; + + sig11 = 1./fu1.^2 + 1./fu2.^2; + sig12 = -1./fu1.^2 + 1./fu2.^2; + sigx = sig11 - sig12.^2./sig11; + + w1p = ntgz - sig12./sig11.*ntgu; + if (largescale) + h11pfun = @(z) sigx.*z - (1/fe)*At(A(z)) + 1/fe^2*(atr'*z)*atr; + [dx, cgres, cgiter] = cgsolve(h11pfun, w1p, cgtol, cgmaxiter, 0); + if (cgres > 1/2) + disp('Cannot solve system. Returning previous iterate. (See Section 4 of notes for more information.)'); + xp = x; up = u; + return + end + Adx = A(dx); + else + H11p = diag(sigx) - (1/fe)*AtA + (1/fe)^2*atr*atr'; + opts.POSDEF = true; opts.SYM = true; + [dx,hcond] = linsolve(H11p, w1p, opts); + if (hcond < 1e-14) + disp('Matrix ill-conditioned. Returning previous iterate. (See Section 4 of notes for more information.)'); + xp = x; up = u; + return + end + Adx = A*dx; + end + du = (1./sig11).*ntgu - (sig12./sig11).*dx; + + % minimum step size that stays in the interior + ifu1 = find((dx-du) > 0); ifu2 = find((-dx-du) > 0); + aqe = Adx'*Adx; bqe = 2*r'*Adx; cqe = r'*r - epsilon^2; + smax = min(1,min([... + -fu1(ifu1)./(dx(ifu1)-du(ifu1)); -fu2(ifu2)./(-dx(ifu2)-du(ifu2)); ... + (-bqe+sqrt(bqe^2-4*aqe*cqe))/(2*aqe) + ])); + s = (0.99)*smax; + + % backtracking line search + suffdec = 0; + backiter = 0; + while (~suffdec) + xp = x + s*dx; up = u + s*du; rp = r + s*Adx; + fu1p = xp - up; fu2p = -xp - up; fep = 1/2*(rp'*rp - epsilon^2); + fp = sum(up) - (1/tau)*(sum(log(-fu1p)) + sum(log(-fu2p)) + log(-fep)); + flin = f + alpha*s*(gradf'*[dx; du]); + suffdec = (fp <= flin); + s = beta*s; + backiter = backiter + 1; + if (backiter > 32) + disp('Stuck on backtracking line search, returning previous iterate. (See Section 4 of notes for more information.)'); + xp = x; up = u; + return + end + end + + % set up for next iteration + x = xp; u = up; r = rp; + fu1 = fu1p; fu2 = fu2p; fe = fep; f = fp; + + lambda2 = -(gradf'*[dx; du]); + stepsize = s*norm([dx; du]); + niter = niter + 1; + done = (lambda2/2 < newtontol) | (niter >= newtonmaxiter); + + disp(sprintf('Newton iter = %d, Functional = %8.3f, Newton decrement = %8.3f, Stepsize = %8.3e', ... + niter, f, lambda2/2, stepsize)); + if (largescale) + disp(sprintf(' CG Res = %8.3e, CG Iter = %d', cgres, cgiter)); + else + disp(sprintf(' H11p condition number = %8.3e', hcond)); + end + +end + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/OMP/greed_omp_qr.m Thu Oct 20 19:36:24 2011 +0000 @@ -0,0 +1,472 @@ +function [s, err_mse, iter_time]=greed_omp_qr(x,A,m,varargin) +% greed_omp_qr: Orthogonal Matching Pursuit algorithm based on QR +% factorisation +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Usage +% [s, err_mse, iter_time]=greed_omp_qr(x,P,m,'option_name','option_value') +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Input +% Mandatory: +% x Observation vector to be decomposed +% P Either: +% 1) An nxm matrix (n must be dimension of x) +% 2) A function handle (type "help function_format" +% for more information) +% Also requires specification of P_trans option. +% 3) An object handle (type "help object_format" for +% more information) +% m length of s +% +% Possible additional options: +% (specify as many as you want using 'option_name','option_value' pairs) +% See below for explanation of options: +%__________________________________________________________________________ +% option_name | available option_values | default +%-------------------------------------------------------------------------- +% stopCrit | M, corr, mse, mse_change | M +% stopTol | number (see below) | n/4 +% P_trans | function_handle (see below) | +% maxIter | positive integer (see below) | n +% verbose | true, false | false +% start_val | vector of length m | zeros +% +% Available stopping criteria : +% M - Extracts exactly M = stopTol elements. +% corr - Stops when maximum correlation between +% residual and atoms is below stopTol value. +% mse - Stops when mean squared error of residual +% is below stopTol value. +% mse_change - Stops when the change in the mean squared +% error falls below stopTol value. +% +% stopTol: Value for stopping criterion. +% +% P_trans: If P is a function handle, then P_trans has to be specified and +% must be a function handle. +% +% maxIter: Maximum number of allowed iterations. +% +% verbose: Logical value to allow algorithm progress to be displayed. +% +% start_val: Allows algorithms to start from partial solution. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Outputs +% s Solution vector +% err_mse Vector containing mse of approximation error for each +% iteration +% iter_time Vector containing computation times for each iteration +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Description +% greed_omp_qr performs a greedy signal decomposition. +% In each iteration a new element is selected depending on the inner +% product between the current residual and columns in P. +% The non-zero elements of s are approximated by orthogonally projecting +% x onto the selected elements in each iteration. +% The algorithm uses QR decomposition. +% +% See Also +% greed_omp_chol, greed_omp_cg, greed_omp_cgp, greed_omp_pinv, +% greed_omp_linsolve, greed_gp, greed_nomp +% +% Copyright (c) 2007 Thomas Blumensath +% +% The University of Edinburgh +% Email: thomas.blumensath@ed.ac.uk +% Comments and bug reports welcome +% +% This file is part of sparsity Version 0.1 +% Created: April 2007 +% +% Part of this toolbox was developed with the support of EPSRC Grant +% D000246/1 +% +% Please read COPYRIGHT.m for terms and conditions. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Default values and initialisation +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +[n1 n2]=size(x); +if n2 == 1 + n=n1; +elseif n1 == 1 + x=x'; + n=n2; +else + display('x must be a vector.'); + return +end + +sigsize = x'*x/n; +initial_given=0; +err_mse = []; +iter_time = []; +STOPCRIT = 'M'; +STOPTOL = ceil(n/4); +MAXITER = n; +verbose = false; +s_initial = zeros(m,1); + + +if verbose + display('Initialising...') +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Output variables +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +switch nargout + case 3 + comp_err=true; + comp_time=true; + case 2 + comp_err=true; + comp_time=false; + case 1 + comp_err=false; + comp_time=false; + case 0 + error('Please assign output variable.') + otherwise + error('Too many output arguments specified') +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Look through options +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% Put option into nice format +Options={}; +OS=nargin-3; +c=1; +for i=1:OS + if isa(varargin{i},'cell') + CellSize=length(varargin{i}); + ThisCell=varargin{i}; + for j=1:CellSize + Options{c}=ThisCell{j}; + c=c+1; + end + else + Options{c}=varargin{i}; + c=c+1; + end +end +OS=length(Options); +if rem(OS,2) + error('Something is wrong with argument name and argument value pairs.') +end + +for i=1:2:OS + switch Options{i} + case {'stopCrit'} + if (strmatch(Options{i+1},{'M'; 'corr'; 'mse'; 'mse_change'},'exact')); + STOPCRIT = Options{i+1}; + else error('stopCrit must be char string [M, corr, mse, mse_change]. Exiting.'); end + case {'stopTol'} + if isa(Options{i+1},'numeric') ; STOPTOL = Options{i+1}; + else error('stopTol must be number. Exiting.'); end + case {'P_trans'} + if isa(Options{i+1},'function_handle'); Pt = Options{i+1}; + else error('P_trans must be function _handle. Exiting.'); end + case {'maxIter'} + if isa(Options{i+1},'numeric'); MAXITER = Options{i+1}; + else error('maxIter must be a number. Exiting.'); end + case {'verbose'} + if isa(Options{i+1},'logical'); verbose = Options{i+1}; + else error('verbose must be a logical. Exiting.'); end + case {'start_val'} + if isa(Options{i+1},'numeric') & length(Options{i+1}) == m ; + s_initial = Options{i+1}; + initial_given=1; + else error('start_val must be a vector of length m. Exiting.'); end + otherwise + error('Unrecognised option. Exiting.') + end +end + + + +if strcmp(STOPCRIT,'M') + maxM=STOPTOL; +else + maxM=MAXITER; +end + +if nargout >=2 + err_mse = zeros(maxM,1); +end +if nargout ==3 + iter_time = zeros(maxM,1); +end + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Make P and Pt functions +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +if isa(A,'float') P =@(z) A*z; Pt =@(z) A'*z; +elseif isobject(A) P =@(z) A*z; Pt =@(z) A'*z; +elseif isa(A,'function_handle') + try + if isa(Pt,'function_handle'); P=A; + else error('If P is a function handle, Pt also needs to be a function handle. Exiting.'); end + catch error('If P is a function handle, Pt needs to be specified. Exiting.'); end +else error('P is of unsupported type. Use matrix, function_handle or object. Exiting.'); end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Random Check to see if dictionary is normalised +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % Commented by Nic, 13 Sept 2011 + % Don't want any slow text output +% mask=zeros(m,1); +% mask(ceil(rand*m))=1; +% nP=norm(P(mask)); +% if abs(1-nP)>1e-3; +% display('Dictionary appears not to have unit norm columns.') +% end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Check if we have enough memory and initialise +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + try Q=zeros(n,maxM); + catch error('Variable size is too large. Please try greed_omp_chol algorithm or reduce MAXITER.') + end + try R=zeros(maxM); + catch error('Variable size is too large. Please try greed_omp_chol algorithm or reduce MAXITER.') + end + + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Do we start from zero or not? +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +if initial_given ==1; + IN = find(s_initial); + if ~isempty(IN) + Residual = x-P(s_initial); + lengthIN=length(IN); + z=[]; + for k=1:length(IN) + % Extract new element + mask=zeros(m,1); + mask(IN(k))=1; + new_element=P(mask); + + % Orthogonalise new element + qP=Q(:,1:k-1)'*new_element; + q=new_element-Q(:,1:k-1)*(qP); + + nq=norm(q); + q=q/nq; + % Update QR factorisation + R(1:k-1,k)=qP; + R(k,k)=nq; + Q(:,k)=q; + + z(k)=q'*x; + end + s = s_initial; + Residual=x-Q(:,k)*z; + oldERR = Residual'*Residual/n; + else + IN = []; + Residual = x; + s = s_initial; + sigsize = x'*x/n; + oldERR = sigsize; + k=0; + z=[]; + end + +else + IN = []; + Residual = x; + s = s_initial; + sigsize = x'*x/n; + oldERR = sigsize; + k=0; + z=[]; +end + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Main algorithm +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +if verbose + display('Main iterations...') +end +tic +t=0; +DR=Pt(Residual); +done = 0; +iter=1; + +%Nic: find norms of dictionary atoms, for RandOMP +% for i = 1:m +% norms(i) = norm(P([zeros(i-1,1) ; 1 ; zeros(m-i,1)])); +% end + +while ~done + + % Select new element + DR(IN)=0; + % Nic: replace selection with random variable + % i.e. Randomized OMP!! + % DON'T FORGET ABOUT THIS!! + [v I]=max(abs(DR)); + %I = randp(exp(abs(DR).^2 ./ (norms.^2)'), [1 1]); + IN=[IN I]; + + + k=k+1; + % Extract new element + mask=zeros(m,1); + mask(IN(k))=1; + new_element=P(mask); + + % Orthogonalise new element + qP=Q(:,1:k-1)'*new_element; + q=new_element-Q(:,1:k-1)*(qP); + + nq=norm(q); + q=q/nq; + % Update QR factorisation + R(1:k-1,k)=qP; + R(k,k)=nq; + Q(:,k)=q; + + z(k)=q'*x; + + % New residual + Residual=Residual-q*(z(k)); + DR=Pt(Residual); + + ERR=Residual'*Residual/n; + if comp_err + err_mse(iter)=ERR; + end + + if comp_time + iter_time(iter)=toc; + end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Are we done yet? +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + if strcmp(STOPCRIT,'M') + if iter >= STOPTOL + done =1; + elseif verbose && toc-t>10 + display(sprintf('Iteration %i. --- %i iterations to go',iter ,STOPTOL-iter)) + t=toc; + end + elseif strcmp(STOPCRIT,'mse') + if comp_err + if err_mse(iter)<STOPTOL; + done = 1; + elseif verbose && toc-t>10 + display(sprintf('Iteration %i. --- %i mse',iter ,err_mse(iter))) + t=toc; + end + else + if ERR<STOPTOL; + done = 1; + elseif verbose && toc-t>10 + display(sprintf('Iteration %i. --- %i mse',iter ,ERR)) + t=toc; + end + end + elseif strcmp(STOPCRIT,'mse_change') && iter >=2 + if comp_err && iter >=2 + if ((err_mse(iter-1)-err_mse(iter))/sigsize <STOPTOL); + done = 1; + elseif verbose && toc-t>10 + display(sprintf('Iteration %i. --- %i mse change',iter ,(err_mse(iter-1)-err_mse(iter))/sigsize )) + t=toc; + end + else + if ((oldERR - ERR)/sigsize < STOPTOL); + done = 1; + elseif verbose && toc-t>10 + display(sprintf('Iteration %i. --- %i mse change',iter ,(oldERR - ERR)/sigsize)) + t=toc; + end + end + elseif strcmp(STOPCRIT,'corr') + if max(abs(DR)) < STOPTOL; + done = 1; + elseif verbose && toc-t>10 + display(sprintf('Iteration %i. --- %i corr',iter ,max(abs(DR)))) + t=toc; + end + end + + % Also stop if residual gets too small or maxIter reached + if comp_err + if err_mse(iter)<1e-16 + display('Stopping. Exact signal representation found!') + done=1; + end + else + + + if iter>1 + if ERR<1e-16 + display('Stopping. Exact signal representation found!') + done=1; + end + end + end + + if iter >= MAXITER + display('Stopping. Maximum number of iterations reached!') + done = 1; + end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% If not done, take another round +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + if ~done + iter=iter+1; + oldERR=ERR; + end +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Now we can solve for s by back-substitution +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + s(IN)=R(1:k,1:k)\z(1:k)'; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Only return as many elements as iterations +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +if nargout >=2 + err_mse = err_mse(1:iter); +end +if nargout ==3 + iter_time = iter_time(1:iter); +end +if verbose + display('Done') +end + +% Change history +% +% 8 of Februray: Algo does no longer stop if dictionary is not normaliesd.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/OMP/omp.py Thu Oct 20 19:36:24 2011 +0000 @@ -0,0 +1,557 @@ +"""Orthogonal matching pursuit algorithms +""" + +# Author: Vlad Niculae +# +# License: BSD Style. + +from warnings import warn + +import numpy as np +from scipy import linalg +from scipy.linalg.lapack import get_lapack_funcs + +#from .base import LinearModel +#from ..utils.arrayfuncs import solve_triangular +from sklearn.linear_model.base import LinearModel +from sklearn.utils.arrayfuncs import solve_triangular + +premature = """ Orthogonal matching pursuit ended prematurely due to linear +dependence in the dictionary. The requested precision might not have been met. +""" + + +def _cholesky_omp(X, y, n_nonzero_coefs, tol=None, copy_X=True): + """Orthogonal Matching Pursuit step using the Cholesky decomposition. + +Parameters: +----------- +X: array, shape = (n_samples, n_features) +Input dictionary. Columns are assumed to have unit norm. + +y: array, shape = (n_samples,) +Input targets + +n_nonzero_coefs: int +Targeted number of non-zero elements + +tol: float +Targeted squared error, if not None overrides n_nonzero_coefs. + +copy_X: bool, optional +Whether the design matrix X must be copied by the algorithm. A false +value is only helpful if X is already Fortran-ordered, otherwise a +copy is made anyway. + +Returns: +-------- +gamma: array, shape = (n_nonzero_coefs,) +Non-zero elements of the solution + +idx: array, shape = (n_nonzero_coefs,) +Indices of the positions of the elements in gamma within the solution +vector + +""" + if copy_X: + X = X.copy('F') + else: # even if we are allowed to overwrite, still copy it if bad order + X = np.asfortranarray(X) + + min_float = np.finfo(X.dtype).eps + nrm2, swap = linalg.get_blas_funcs(('nrm2', 'swap'), (X,)) + potrs, = get_lapack_funcs(('potrs',), (X,)) + + alpha = np.dot(X.T, y) + residual = y + n_active = 0 + indices = range(X.shape[1]) # keeping track of swapping + + max_features = X.shape[1] if tol is not None else n_nonzero_coefs + L = np.empty((max_features, max_features), dtype=X.dtype) + L[0, 0] = 1. + + while True: + lam = np.argmax(np.abs(np.dot(X.T, residual))) + if lam < n_active or alpha[lam] ** 2 < min_float: + # atom already selected or inner product too small + warn(premature) + break + if n_active > 0: + # Updates the Cholesky decomposition of X' X + L[n_active, :n_active] = np.dot(X[:, :n_active].T, X[:, lam]) + solve_triangular(L[:n_active, :n_active], L[n_active, :n_active]) + v = nrm2(L[n_active, :n_active]) ** 2 + if 1 - v <= min_float: # selected atoms are dependent + warn(premature) + break + L[n_active, n_active] = np.sqrt(1 - v) + X.T[n_active], X.T[lam] = swap(X.T[n_active], X.T[lam]) + alpha[n_active], alpha[lam] = alpha[lam], alpha[n_active] + indices[n_active], indices[lam] = indices[lam], indices[n_active] + n_active += 1 + # solves LL'x = y as a composition of two triangular systems + gamma, _ = potrs(L[:n_active, :n_active], alpha[:n_active], lower=True, + overwrite_b=False) + + residual = y - np.dot(X[:, :n_active], gamma) + if tol is not None and nrm2(residual) ** 2 <= tol: + break + elif n_active == max_features: + break + + return gamma, indices[:n_active] + + +def _gram_omp(Gram, Xy, n_nonzero_coefs, tol_0=None, tol=None, + copy_Gram=True, copy_Xy=True): + """Orthogonal Matching Pursuit step on a precomputed Gram matrix. + +This function uses the the Cholesky decomposition method. + +Parameters: +----------- +Gram: array, shape = (n_features, n_features) +Gram matrix of the input data matrix + +Xy: array, shape = (n_features,) +Input targets + +n_nonzero_coefs: int +Targeted number of non-zero elements + +tol_0: float +Squared norm of y, required if tol is not None. + +tol: float +Targeted squared error, if not None overrides n_nonzero_coefs. + +copy_Gram: bool, optional +Whether the gram matrix must be copied by the algorithm. A false +value is only helpful if it is already Fortran-ordered, otherwise a +copy is made anyway. + +copy_Xy: bool, optional +Whether the covariance vector Xy must be copied by the algorithm. +If False, it may be overwritten. + +Returns: +-------- +gamma: array, shape = (n_nonzero_coefs,) +Non-zero elements of the solution + +idx: array, shape = (n_nonzero_coefs,) +Indices of the positions of the elements in gamma within the solution +vector + +""" + Gram = Gram.copy('F') if copy_Gram else np.asfortranarray(Gram) + + if copy_Xy: + Xy = Xy.copy() + + min_float = np.finfo(Gram.dtype).eps + nrm2, swap = linalg.get_blas_funcs(('nrm2', 'swap'), (Gram,)) + potrs, = get_lapack_funcs(('potrs',), (Gram,)) + + indices = range(len(Gram)) # keeping track of swapping + alpha = Xy + tol_curr = tol_0 + delta = 0 + n_active = 0 + + max_features = len(Gram) if tol is not None else n_nonzero_coefs + L = np.empty((max_features, max_features), dtype=Gram.dtype) + L[0, 0] = 1. + + while True: + lam = np.argmax(np.abs(alpha)) + if lam < n_active or alpha[lam] ** 2 < min_float: + # selected same atom twice, or inner product too small + warn(premature) + break + if n_active > 0: + L[n_active, :n_active] = Gram[lam, :n_active] + solve_triangular(L[:n_active, :n_active], L[n_active, :n_active]) + v = nrm2(L[n_active, :n_active]) ** 2 + if 1 - v <= min_float: # selected atoms are dependent + warn(premature) + break + L[n_active, n_active] = np.sqrt(1 - v) + Gram[n_active], Gram[lam] = swap(Gram[n_active], Gram[lam]) + Gram.T[n_active], Gram.T[lam] = swap(Gram.T[n_active], Gram.T[lam]) + indices[n_active], indices[lam] = indices[lam], indices[n_active] + Xy[n_active], Xy[lam] = Xy[lam], Xy[n_active] + n_active += 1 + # solves LL'x = y as a composition of two triangular systems + gamma, _ = potrs(L[:n_active, :n_active], Xy[:n_active], lower=True, + overwrite_b=False) + + beta = np.dot(Gram[:, :n_active], gamma) + alpha = Xy - beta + if tol is not None: + tol_curr += delta + delta = np.inner(gamma, beta[:n_active]) + tol_curr -= delta + if tol_curr <= tol: + break + elif n_active == max_features: + break + + return gamma, indices[:n_active] + + +def orthogonal_mp(X, y, n_nonzero_coefs=None, tol=None, precompute_gram=False, + copy_X=True): + """Orthogonal Matching Pursuit (OMP) + +Solves n_targets Orthogonal Matching Pursuit problems. +An instance of the problem has the form: + +When parametrized by the number of non-zero coefficients using +`n_nonzero_coefs`: +argmin ||y - X\gamma||^2 subject to ||\gamma||_0 <= n_{nonzero coefs} + +When parametrized by error using the parameter `tol`: +argmin ||\gamma||_0 subject to ||y - X\gamma||^2 <= tol + +Parameters +---------- +X: array, shape = (n_samples, n_features) +Input data. Columns are assumed to have unit norm. + +y: array, shape = (n_samples,) or (n_samples, n_targets) +Input targets + +n_nonzero_coefs: int +Desired number of non-zero entries in the solution. If None (by +default) this value is set to 10% of n_features. + +tol: float +Maximum norm of the residual. If not None, overrides n_nonzero_coefs. + +precompute_gram: {True, False, 'auto'}, +Whether to perform precomputations. Improves performance when n_targets +or n_samples is very large. + +copy_X: bool, optional +Whether the design matrix X must be copied by the algorithm. A false +value is only helpful if X is already Fortran-ordered, otherwise a +copy is made anyway. + +Returns +------- +coef: array, shape = (n_features,) or (n_features, n_targets) +Coefficients of the OMP solution + +See also +-------- +OrthogonalMatchingPursuit +orthogonal_mp_gram +lars_path +decomposition.sparse_encode +decomposition.sparse_encode_parallel + +Notes +----- +Orthogonal matching pursuit was introduced in G. Mallat, Z. Zhang, +Matching pursuits with time-frequency dictionaries, IEEE Transactions on +Signal Processing, Vol. 41, No. 12. (December 1993), pp. 3397-3415. +(http://blanche.polytechnique.fr/~mallat/papiers/MallatPursuit93.pdf) + +This implementation is based on Rubinstein, R., Zibulevsky, M. and Elad, +M., Efficient Implementation of the K-SVD Algorithm using Batch Orthogonal +Matching Pursuit Technical Report - CS Technion, April 2008. +http://www.cs.technion.ac.il/~ronrubin/Publications/KSVX-OMP-v2.pdf + +""" + X, y = map(np.asanyarray, (X, y)) + if y.ndim == 1: + y = y[:, np.newaxis] + if copy_X: + X = X.copy('F') + copy_X = False + else: + X = np.asfortranarray(X) + if y.shape[1] > 1: # subsequent targets will be affected + copy_X = True + if n_nonzero_coefs == None and tol == None: + n_nonzero_coefs = int(0.1 * X.shape[1]) + if tol is not None and tol < 0: + raise ValueError("Epsilon cannot be negative") + if tol is None and n_nonzero_coefs <= 0: + raise ValueError("The number of atoms must be positive") + if tol is None and n_nonzero_coefs > X.shape[1]: + raise ValueError("The number of atoms cannot be more than the number \ +of features") + if precompute_gram == 'auto': + precompute_gram = X.shape[0] > X.shape[1] + if precompute_gram: + G = np.dot(X.T, X) + G = np.asfortranarray(G) + Xy = np.dot(X.T, y) + if tol is not None: + norms_squared = np.sum((y ** 2), axis=0) + else: + norms_squared = None + return orthogonal_mp_gram(G, Xy, n_nonzero_coefs, tol, norms_squared, + copy_Gram=copy_X, copy_Xy=False) + + coef = np.zeros((X.shape[1], y.shape[1])) + for k in xrange(y.shape[1]): + x, idx = _cholesky_omp(X, y[:, k], n_nonzero_coefs, tol, + copy_X=copy_X) + coef[idx, k] = x + return np.squeeze(coef) + + +def orthogonal_mp_gram(Gram, Xy, n_nonzero_coefs=None, tol=None, + norms_squared=None, copy_Gram=True, + copy_Xy=True): + """Gram Orthogonal Matching Pursuit (OMP) + +Solves n_targets Orthogonal Matching Pursuit problems using only +the Gram matrix X.T * X and the product X.T * y. + +Parameters +---------- +Gram: array, shape = (n_features, n_features) +Gram matrix of the input data: X.T * X + +Xy: array, shape = (n_features,) or (n_features, n_targets) +Input targets multiplied by X: X.T * y + +n_nonzero_coefs: int +Desired number of non-zero entries in the solution. If None (by +default) this value is set to 10% of n_features. + +tol: float +Maximum norm of the residual. If not None, overrides n_nonzero_coefs. + +norms_squared: array-like, shape = (n_targets,) +Squared L2 norms of the lines of y. Required if tol is not None. + +copy_Gram: bool, optional +Whether the gram matrix must be copied by the algorithm. A false +value is only helpful if it is already Fortran-ordered, otherwise a +copy is made anyway. + +copy_Xy: bool, optional +Whether the covariance vector Xy must be copied by the algorithm. +If False, it may be overwritten. + +Returns +------- +coef: array, shape = (n_features,) or (n_features, n_targets) +Coefficients of the OMP solution + +See also +-------- +OrthogonalMatchingPursuit +orthogonal_mp +lars_path +decomposition.sparse_encode +decomposition.sparse_encode_parallel + +Notes +----- +Orthogonal matching pursuit was introduced in G. Mallat, Z. Zhang, +Matching pursuits with time-frequency dictionaries, IEEE Transactions on +Signal Processing, Vol. 41, No. 12. (December 1993), pp. 3397-3415. +(http://blanche.polytechnique.fr/~mallat/papiers/MallatPursuit93.pdf) + +This implementation is based on Rubinstein, R., Zibulevsky, M. and Elad, +M., Efficient Implementation of the K-SVD Algorithm using Batch Orthogonal +Matching Pursuit Technical Report - CS Technion, April 2008. +http://www.cs.technion.ac.il/~ronrubin/Publications/KSVX-OMP-v2.pdf + +""" + Gram, Xy = map(np.asanyarray, (Gram, Xy)) + if Xy.ndim == 1: + Xy = Xy[:, np.newaxis] + if tol is not None: + norms_squared = [norms_squared] + + if n_nonzero_coefs == None and tol is None: + n_nonzero_coefs = int(0.1 * len(Gram)) + if tol is not None and norms_squared == None: + raise ValueError('Gram OMP needs the precomputed norms in order \ +to evaluate the error sum of squares.') + if tol is not None and tol < 0: + raise ValueError("Epsilon cennot be negative") + if tol is None and n_nonzero_coefs <= 0: + raise ValueError("The number of atoms must be positive") + if tol is None and n_nonzero_coefs > len(Gram): + raise ValueError("The number of atoms cannot be more than the number \ +of features") + coef = np.zeros((len(Gram), Xy.shape[1])) + for k in range(Xy.shape[1]): + x, idx = _gram_omp(Gram, Xy[:, k], n_nonzero_coefs, + norms_squared[k] if tol is not None else None, tol, + copy_Gram=copy_Gram, copy_Xy=copy_Xy) + coef[idx, k] = x + return np.squeeze(coef) + + +class OrthogonalMatchingPursuit(LinearModel): + """Orthogonal Mathching Pursuit model (OMP) + +Parameters +---------- +n_nonzero_coefs: int, optional +Desired number of non-zero entries in the solution. If None (by +default) this value is set to 10% of n_features. + +tol: float, optional +Maximum norm of the residual. If not None, overrides n_nonzero_coefs. + +fit_intercept: boolean, optional +whether to calculate the intercept for this model. If set +to false, no intercept will be used in calculations +(e.g. data is expected to be already centered). + +normalize: boolean, optional +If False, the regressors X are assumed to be already normalized. + +precompute_gram: {True, False, 'auto'}, +Whether to use a precomputed Gram and Xy matrix to speed up +calculations. Improves performance when `n_targets` or `n_samples` is +very large. Note that if you already have such matrices, you can pass +them directly to the fit method. + +copy_X: bool, optional +Whether the design matrix X must be copied by the algorithm. A false +value is only helpful if X is already Fortran-ordered, otherwise a +copy is made anyway. + +copy_Gram: bool, optional +Whether the gram matrix must be copied by the algorithm. A false +value is only helpful if X is already Fortran-ordered, otherwise a +copy is made anyway. + +copy_Xy: bool, optional +Whether the covariance vector Xy must be copied by the algorithm. +If False, it may be overwritten. + + +Attributes +---------- +coef_: array, shape = (n_features,) or (n_features, n_targets) +parameter vector (w in the fomulation formula) + +intercept_: float or array, shape =(n_targets,) +independent term in decision function. + +Notes +----- +Orthogonal matching pursuit was introduced in G. Mallat, Z. Zhang, +Matching pursuits with time-frequency dictionaries, IEEE Transactions on +Signal Processing, Vol. 41, No. 12. (December 1993), pp. 3397-3415. +(http://blanche.polytechnique.fr/~mallat/papiers/MallatPursuit93.pdf) + +This implementation is based on Rubinstein, R., Zibulevsky, M. and Elad, +M., Efficient Implementation of the K-SVD Algorithm using Batch Orthogonal +Matching Pursuit Technical Report - CS Technion, April 2008. +http://www.cs.technion.ac.il/~ronrubin/Publications/KSVX-OMP-v2.pdf + +See also +-------- +orthogonal_mp +orthogonal_mp_gram +lars_path +Lars +LassoLars +decomposition.sparse_encode +decomposition.sparse_encode_parallel + +""" + def __init__(self, copy_X=True, copy_Gram=True, + copy_Xy=True, n_nonzero_coefs=None, tol=None, + fit_intercept=True, normalize=True, precompute_gram=False): + self.n_nonzero_coefs = n_nonzero_coefs + self.tol = tol + self.fit_intercept = fit_intercept + self.normalize = normalize + self.precompute_gram = precompute_gram + self.copy_Gram = copy_Gram + self.copy_Xy = copy_Xy + self.copy_X = copy_X + + def fit(self, X, y, Gram=None, Xy=None): + """Fit the model using X, y as training data. + +Parameters +---------- +X: array-like, shape = (n_samples, n_features) +Training data. + +y: array-like, shape = (n_samples,) or (n_samples, n_targets) +Target values. + +Gram: array-like, shape = (n_features, n_features) (optional) +Gram matrix of the input data: X.T * X + +Xy: array-like, shape = (n_features,) or (n_features, n_targets) +(optional) +Input targets multiplied by X: X.T * y + + +Returns +------- +self: object +returns an instance of self. +""" + X = np.atleast_2d(X) + y = np.atleast_1d(y) + n_features = X.shape[1] + + X, y, X_mean, y_mean, X_std = self._center_data(X, y, + self.fit_intercept, + self.normalize, + self.copy_X) + + if self.n_nonzero_coefs == None and self.tol is None: + self.n_nonzero_coefs = int(0.1 * n_features) + + if Gram is not None: + Gram = np.atleast_2d(Gram) + + if self.copy_Gram: + copy_Gram = False + Gram = Gram.copy('F') + else: + Gram = np.asfortranarray(Gram) + + copy_Gram = self.copy_Gram + if y.shape[1] > 1: # subsequent targets will be affected + copy_Gram = True + + if Xy is None: + Xy = np.dot(X.T, y) + else: + if self.copy_Xy: + Xy = Xy.copy() + if self.normalize: + if len(Xy.shape) == 1: + Xy /= X_std + else: + Xy /= X_std[:, np.newaxis] + + if self.normalize: + Gram /= X_std + Gram /= X_std[:, np.newaxis] + + norms_sq = np.sum(y ** 2, axis=0) if self.tol is not None else None + self.coef_ = orthogonal_mp_gram(Gram, Xy, self.n_nonzero_coefs, + self.tol, norms_sq, + copy_Gram, True).T + else: + precompute_gram = self.precompute_gram + if precompute_gram == 'auto': + precompute_gram = X.shape[0] > X.shape[1] + self.coef_ = orthogonal_mp(X, y, self.n_nonzero_coefs, self.tol, + precompute_gram=self.precompute_gram, + copy_X=self.copy_X).T + + self._set_intercept(X_mean, y_mean, X_std) + return self \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/OMP/omp_downloaded.py Thu Oct 20 19:36:24 2011 +0000 @@ -0,0 +1,558 @@ +"""Orthogonal matching pursuit algorithms +""" + +# Author: Vlad Niculae +# +# License: BSD Style. + +from warnings import warn + +import numpy as np +from scipy import linalg +from scipy.linalg.lapack import get_lapack_funcs + +from sklearn.linear_model.base import LinearModel +from sklearn.utils.arrayfuncs import solve_triangular + +premature = """ Orthogonal matching pursuit ended prematurely due to linear +dependence in the dictionary. The requested precision might not have been met. +""" + + +def _cholesky_omp(X, y, n_nonzero_coefs, tol=None, copy_X=True): + """Orthogonal Matching Pursuit step using the Cholesky decomposition. + + Parameters: + ----------- + X: array, shape = (n_samples, n_features) + Input dictionary. Columns are assumed to have unit norm. + + y: array, shape = (n_samples,) + Input targets + + n_nonzero_coefs: int + Targeted number of non-zero elements + + tol: float + Targeted squared error, if not None overrides n_nonzero_coefs. + + copy_X: bool, optional + Whether the design matrix X must be copied by the algorithm. A false + value is only helpful if X is already Fortran-ordered, otherwise a + copy is made anyway. + + Returns: + -------- + gamma: array, shape = (n_nonzero_coefs,) + Non-zero elements of the solution + + idx: array, shape = (n_nonzero_coefs,) + Indices of the positions of the elements in gamma within the solution + vector + + """ + if copy_X: + X = X.copy('F') + else: # even if we are allowed to overwrite, still copy it if bad order + X = np.asfortranarray(X) + + min_float = np.finfo(X.dtype).eps + nrm2, swap = linalg.get_blas_funcs(('nrm2', 'swap'), (X,)) + potrs, = get_lapack_funcs(('potrs',), (X,)) + + alpha = np.dot(X.T, y) + residual = y + n_active = 0 + idx = [] + + max_features = X.shape[1] if tol is not None else n_nonzero_coefs + L = np.empty((max_features, max_features), dtype=X.dtype) + L[0, 0] = 1. + + while True: + lam = np.argmax(np.abs(np.dot(X.T, residual))) + if lam < n_active or alpha[lam] ** 2 < min_float: + # atom already selected or inner product too small + warn(premature) + break + if n_active > 0: + # Updates the Cholesky decomposition of X' X + L[n_active, :n_active] = np.dot(X[:, :n_active].T, X[:, lam]) + solve_triangular(L[:n_active, :n_active], L[n_active, :n_active]) + v = nrm2(L[n_active, :n_active]) ** 2 + if 1 - v <= min_float: # selected atoms are dependent + warn(premature) + break + L[n_active, n_active] = np.sqrt(1 - v) + idx.append(lam) + X.T[n_active], X.T[lam] = swap(X.T[n_active], X.T[lam]) + alpha[n_active], alpha[lam] = alpha[lam], alpha[n_active] + n_active += 1 + # solves LL'x = y as a composition of two triangular systems + gamma, _ = potrs(L[:n_active, :n_active], alpha[:n_active], lower=True, + overwrite_b=False) + + residual = y - np.dot(X[:, :n_active], gamma) + if tol is not None and nrm2(residual) ** 2 <= tol: + break + elif n_active == max_features: + break + + return gamma, idx + + +def _gram_omp(Gram, Xy, n_nonzero_coefs, tol_0=None, tol=None, + copy_Gram=True, copy_Xy=True): + """Orthogonal Matching Pursuit step on a precomputed Gram matrix. + + This function uses the the Cholesky decomposition method. + + Parameters: + ----------- + Gram: array, shape = (n_features, n_features) + Gram matrix of the input data matrix + + Xy: array, shape = (n_features,) + Input targets + + n_nonzero_coefs: int + Targeted number of non-zero elements + + tol_0: float + Squared norm of y, required if tol is not None. + + tol: float + Targeted squared error, if not None overrides n_nonzero_coefs. + + copy_Gram: bool, optional + Whether the gram matrix must be copied by the algorithm. A false + value is only helpful if it is already Fortran-ordered, otherwise a + copy is made anyway. + + copy_Xy: bool, optional + Whether the covariance vector Xy must be copied by the algorithm. + If False, it may be overwritten. + + Returns: + -------- + gamma: array, shape = (n_nonzero_coefs,) + Non-zero elements of the solution + + idx: array, shape = (n_nonzero_coefs,) + Indices of the positions of the elements in gamma within the solution + vector + + """ + Gram = Gram.copy('F') if copy_Gram else np.asfortranarray(Gram) + + if copy_Xy: + Xy = Xy.copy() + + min_float = np.finfo(Gram.dtype).eps + nrm2, swap = linalg.get_blas_funcs(('nrm2', 'swap'), (Gram,)) + potrs, = get_lapack_funcs(('potrs',), (Gram,)) + + idx = [] + alpha = Xy + tol_curr = tol_0 + delta = 0 + n_active = 0 + + max_features = len(Gram) if tol is not None else n_nonzero_coefs + L = np.empty((max_features, max_features), dtype=Gram.dtype) + L[0, 0] = 1. + + while True: + lam = np.argmax(np.abs(alpha)) + if lam < n_active or alpha[lam] ** 2 < min_float: + # selected same atom twice, or inner product too small + warn(premature) + break + if n_active > 0: + L[n_active, :n_active] = Gram[lam, :n_active] + solve_triangular(L[:n_active, :n_active], L[n_active, :n_active]) + v = nrm2(L[n_active, :n_active]) ** 2 + if 1 - v <= min_float: # selected atoms are dependent + warn(premature) + break + L[n_active, n_active] = np.sqrt(1 - v) + idx.append(lam) + Gram[n_active], Gram[lam] = swap(Gram[n_active], Gram[lam]) + Gram.T[n_active], Gram.T[lam] = swap(Gram.T[n_active], Gram.T[lam]) + Xy[n_active], Xy[lam] = Xy[lam], Xy[n_active] + n_active += 1 + # solves LL'x = y as a composition of two triangular systems + gamma, _ = potrs(L[:n_active, :n_active], Xy[:n_active], lower=True, + overwrite_b=False) + + beta = np.dot(Gram[:, :n_active], gamma) + alpha = Xy - beta + if tol is not None: + tol_curr += delta + delta = np.inner(gamma, beta[:n_active]) + tol_curr -= delta + if tol_curr <= tol: + break + elif n_active == max_features: + break + + return gamma, idx + + +def orthogonal_mp(X, y, n_nonzero_coefs=None, tol=None, precompute_gram=False, + copy_X=True): + """Orthogonal Matching Pursuit (OMP) + + Solves n_targets Orthogonal Matching Pursuit problems. + An instance of the problem has the form: + + When parametrized by the number of non-zero coefficients using + `n_nonzero_coefs`: + argmin ||y - X\gamma||^2 subject to ||\gamma||_0 <= n_{nonzero coefs} + + When parametrized by error using the parameter `tol`: + argmin ||\gamma||_0 subject to ||y - X\gamma||^2 <= tol + + Parameters + ---------- + X: array, shape = (n_samples, n_features) + Input data. Columns are assumed to have unit norm. + + y: array, shape = (n_samples,) or (n_samples, n_targets) + Input targets + + n_nonzero_coefs: int + Desired number of non-zero entries in the solution. If None (by + default) this value is set to 10% of n_features. + + tol: float + Maximum norm of the residual. If not None, overrides n_nonzero_coefs. + + precompute_gram: {True, False, 'auto'}, + Whether to perform precomputations. Improves performance when n_targets + or n_samples is very large. + + copy_X: bool, optional + Whether the design matrix X must be copied by the algorithm. A false + value is only helpful if X is already Fortran-ordered, otherwise a + copy is made anyway. + + Returns + ------- + coef: array, shape = (n_features,) or (n_features, n_targets) + Coefficients of the OMP solution + + See also + -------- + OrthogonalMatchingPursuit + orthogonal_mp_gram + lars_path + decomposition.sparse_encode + decomposition.sparse_encode_parallel + + Notes + ----- + Orthogonal matching pursuit was introduced in G. Mallat, Z. Zhang, + Matching pursuits with time-frequency dictionaries, IEEE Transactions on + Signal Processing, Vol. 41, No. 12. (December 1993), pp. 3397-3415. + (http://blanche.polytechnique.fr/~mallat/papiers/MallatPursuit93.pdf) + + This implementation is based on Rubinstein, R., Zibulevsky, M. and Elad, + M., Efficient Implementation of the K-SVD Algorithm using Batch Orthogonal + Matching Pursuit Technical Report - CS Technion, April 2008. + http://www.cs.technion.ac.il/~ronrubin/Publications/KSVX-OMP-v2.pdf + + """ + + print "New version of OMP" + + X, y = map(np.asanyarray, (X, y)) + if y.ndim == 1: + y = y[:, np.newaxis] + if copy_X: + X = X.copy('F') + copy_X = False + else: + X = np.asfortranarray(X) + if y.shape[1] > 1: # subsequent targets will be affected + copy_X = True + if n_nonzero_coefs == None and tol == None: + n_nonzero_coefs = int(0.1 * X.shape[1]) + if tol is not None and tol < 0: + raise ValueError("Epsilon cannot be negative") + if tol is None and n_nonzero_coefs <= 0: + raise ValueError("The number of atoms must be positive") + if tol is None and n_nonzero_coefs > X.shape[1]: + raise ValueError("The number of atoms cannot be more than the number \ + of features") + if precompute_gram == 'auto': + precompute_gram = X.shape[0] > X.shape[1] + if precompute_gram: + G = np.dot(X.T, X) + G = np.asfortranarray(G) + Xy = np.dot(X.T, y) + if tol is not None: + norms_squared = np.sum((y ** 2), axis=0) + else: + norms_squared = None + return orthogonal_mp_gram(G, Xy, n_nonzero_coefs, tol, norms_squared, + copy_Gram=copy_X, copy_Xy=False) + + coef = np.zeros((X.shape[1], y.shape[1])) + for k in xrange(y.shape[1]): + x, idx = _cholesky_omp(X, y[:, k], n_nonzero_coefs, tol, + copy_X=copy_X) + coef[idx, k] = x + return np.squeeze(coef) + + +def orthogonal_mp_gram(Gram, Xy, n_nonzero_coefs=None, tol=None, + norms_squared=None, copy_Gram=True, + copy_Xy=True): + """Gram Orthogonal Matching Pursuit (OMP) + + Solves n_targets Orthogonal Matching Pursuit problems using only + the Gram matrix X.T * X and the product X.T * y. + + Parameters + ---------- + Gram: array, shape = (n_features, n_features) + Gram matrix of the input data: X.T * X + + Xy: array, shape = (n_features,) or (n_features, n_targets) + Input targets multiplied by X: X.T * y + + n_nonzero_coefs: int + Desired number of non-zero entries in the solution. If None (by + default) this value is set to 10% of n_features. + + tol: float + Maximum norm of the residual. If not None, overrides n_nonzero_coefs. + + norms_squared: array-like, shape = (n_targets,) + Squared L2 norms of the lines of y. Required if tol is not None. + + copy_Gram: bool, optional + Whether the gram matrix must be copied by the algorithm. A false + value is only helpful if it is already Fortran-ordered, otherwise a + copy is made anyway. + + copy_Xy: bool, optional + Whether the covariance vector Xy must be copied by the algorithm. + If False, it may be overwritten. + + Returns + ------- + coef: array, shape = (n_features,) or (n_features, n_targets) + Coefficients of the OMP solution + + See also + -------- + OrthogonalMatchingPursuit + orthogonal_mp + lars_path + decomposition.sparse_encode + decomposition.sparse_encode_parallel + + Notes + ----- + Orthogonal matching pursuit was introduced in G. Mallat, Z. Zhang, + Matching pursuits with time-frequency dictionaries, IEEE Transactions on + Signal Processing, Vol. 41, No. 12. (December 1993), pp. 3397-3415. + (http://blanche.polytechnique.fr/~mallat/papiers/MallatPursuit93.pdf) + + This implementation is based on Rubinstein, R., Zibulevsky, M. and Elad, + M., Efficient Implementation of the K-SVD Algorithm using Batch Orthogonal + Matching Pursuit Technical Report - CS Technion, April 2008. + http://www.cs.technion.ac.il/~ronrubin/Publications/KSVX-OMP-v2.pdf + + """ + Gram, Xy = map(np.asanyarray, (Gram, Xy)) + if Xy.ndim == 1: + Xy = Xy[:, np.newaxis] + if tol is not None: + norms_squared = [norms_squared] + + if n_nonzero_coefs == None and tol is None: + n_nonzero_coefs = int(0.1 * len(Gram)) + if tol is not None and norms_squared == None: + raise ValueError('Gram OMP needs the precomputed norms in order \ + to evaluate the error sum of squares.') + if tol is not None and tol < 0: + raise ValueError("Epsilon cennot be negative") + if tol is None and n_nonzero_coefs <= 0: + raise ValueError("The number of atoms must be positive") + if tol is None and n_nonzero_coefs > len(Gram): + raise ValueError("The number of atoms cannot be more than the number \ + of features") + coef = np.zeros((len(Gram), Xy.shape[1])) + for k in range(Xy.shape[1]): + x, idx = _gram_omp(Gram, Xy[:, k], n_nonzero_coefs, + norms_squared[k] if tol is not None else None, tol, + copy_Gram=copy_Gram, copy_Xy=copy_Xy) + coef[idx, k] = x + return np.squeeze(coef) + + +class OrthogonalMatchingPursuit(LinearModel): + """Orthogonal Mathching Pursuit model (OMP) + + Parameters + ---------- + n_nonzero_coefs: int, optional + Desired number of non-zero entries in the solution. If None (by + default) this value is set to 10% of n_features. + + tol: float, optional + Maximum norm of the residual. If not None, overrides n_nonzero_coefs. + + fit_intercept: boolean, optional + whether to calculate the intercept for this model. If set + to false, no intercept will be used in calculations + (e.g. data is expected to be already centered). + + normalize: boolean, optional + If False, the regressors X are assumed to be already normalized. + + precompute_gram: {True, False, 'auto'}, + Whether to use a precomputed Gram and Xy matrix to speed up + calculations. Improves performance when `n_targets` or `n_samples` is + very large. Note that if you already have such matrices, you can pass + them directly to the fit method. + + copy_X: bool, optional + Whether the design matrix X must be copied by the algorithm. A false + value is only helpful if X is already Fortran-ordered, otherwise a + copy is made anyway. + + copy_Gram: bool, optional + Whether the gram matrix must be copied by the algorithm. A false + value is only helpful if X is already Fortran-ordered, otherwise a + copy is made anyway. + + copy_Xy: bool, optional + Whether the covariance vector Xy must be copied by the algorithm. + If False, it may be overwritten. + + + Attributes + ---------- + coef_: array, shape = (n_features,) or (n_features, n_targets) + parameter vector (w in the fomulation formula) + + intercept_: float or array, shape =(n_targets,) + independent term in decision function. + + Notes + ----- + Orthogonal matching pursuit was introduced in G. Mallat, Z. Zhang, + Matching pursuits with time-frequency dictionaries, IEEE Transactions on + Signal Processing, Vol. 41, No. 12. (December 1993), pp. 3397-3415. + (http://blanche.polytechnique.fr/~mallat/papiers/MallatPursuit93.pdf) + + This implementation is based on Rubinstein, R., Zibulevsky, M. and Elad, + M., Efficient Implementation of the K-SVD Algorithm using Batch Orthogonal + Matching Pursuit Technical Report - CS Technion, April 2008. + http://www.cs.technion.ac.il/~ronrubin/Publications/KSVX-OMP-v2.pdf + + See also + -------- + orthogonal_mp + orthogonal_mp_gram + lars_path + Lars + LassoLars + decomposition.sparse_encode + decomposition.sparse_encode_parallel + + """ + def __init__(self, copy_X=True, copy_Gram=True, + copy_Xy=True, n_nonzero_coefs=None, tol=None, + fit_intercept=True, normalize=True, precompute_gram=False): + self.n_nonzero_coefs = n_nonzero_coefs + self.tol = tol + self.fit_intercept = fit_intercept + self.normalize = normalize + self.precompute_gram = precompute_gram + self.copy_Gram = copy_Gram + self.copy_Xy = copy_Xy + self.copy_X = copy_X + + def fit(self, X, y, Gram=None, Xy=None): + """Fit the model using X, y as training data. + + Parameters + ---------- + X: array-like, shape = (n_samples, n_features) + Training data. + + y: array-like, shape = (n_samples,) or (n_samples, n_targets) + Target values. + + Gram: array-like, shape = (n_features, n_features) (optional) + Gram matrix of the input data: X.T * X + + Xy: array-like, shape = (n_features,) or (n_features, n_targets) + (optional) + Input targets multiplied by X: X.T * y + + + Returns + ------- + self: object + returns an instance of self. + """ + X = np.atleast_2d(X) + y = np.atleast_1d(y) + n_features = X.shape[1] + + X, y, X_mean, y_mean, X_std = self._center_data(X, y, + self.fit_intercept, + self.normalize, + self.copy_X) + + if self.n_nonzero_coefs == None and self.tol is None: + self.n_nonzero_coefs = int(0.1 * n_features) + + if Gram is not None: + Gram = np.atleast_2d(Gram) + + if self.copy_Gram: + copy_Gram = False + Gram = Gram.copy('F') + else: + Gram = np.asfortranarray(Gram) + + copy_Gram = self.copy_Gram + if y.shape[1] > 1: # subsequent targets will be affected + copy_Gram = True + + if Xy is None: + Xy = np.dot(X.T, y) + else: + if self.copy_Xy: + Xy = Xy.copy() + if self.normalize: + if len(Xy.shape) == 1: + Xy /= X_std + else: + Xy /= X_std[:, np.newaxis] + + if self.normalize: + Gram /= X_std + Gram /= X_std[:, np.newaxis] + + norms_sq = np.sum(y ** 2, axis=0) if self.tol is not None else None + self.coef_ = orthogonal_mp_gram(Gram, Xy, self.n_nonzero_coefs, + self.tol, norms_sq, + copy_Gram, True).T + else: + precompute_gram = self.precompute_gram + if precompute_gram == 'auto': + precompute_gram = X.shape[0] > X.shape[1] + self.coef_ = orthogonal_mp(X, y, self.n_nonzero_coefs, self.tol, + precompute_gram=self.precompute_gram, + copy_X=self.copy_X).T + + self._set_intercept(X_mean, y_mean, X_std) + return self
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/OMP/omp_phase3.py Thu Oct 20 19:36:24 2011 +0000 @@ -0,0 +1,986 @@ +""" +#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=# +# Bob L. Sturm <bst@create.aau.dk> 20111018 +# Department of Architecture, Design and Media Technology +# Aalborg University Copenhagen +# Lautrupvang 15, 2750 Ballerup, Denmark +#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=# +""" + +import numpy as np +import scipy.linalg +from sklearn.utils import check_random_state +import time +import math +import omp + +""" +Run a problem suite involving sparse vectors in +ambientDimension dimensional space, with a resolution +in the phase plane of numGradations x numGradations, +and at each indeterminacy and sparsity pair run +numTrials independent trials. + +Outputs a text file denoting successes at each phase point. +For more on phase transitions, see: +D. L. Donoho and J. Tanner, "Precise undersampling theorems," +Proc. IEEE, vol. 98, no. 6, pp. 913-924, June 2010. +""" + +def runProblemSuite(ambientDimension,numGradations,numTrials): + + idx = np.arange(ambientDimension) + phaseDelta = np.linspace(0.05,1,numGradations) + phaseRho = np.linspace(0.05,1,numGradations) + success = np.zeros((numGradations, numGradations)) + + #Nic: init timers + t1all = 0 + t2all = 0 + t3all = 0 + + deltaCounter = 0 + # delta is number of measurements/ + for delta in phaseDelta[:17]: + rhoCounter = 0 + for rho in phaseRho: + print(deltaCounter,rhoCounter) + numMeasurements = int(delta*ambientDimension) + sparsity = int(rho*numMeasurements) + # how do I set the following to be random each time? + generator = check_random_state(100) + # create unit norm dictionary + D = generator.randn(numMeasurements, ambientDimension) + D /= np.sqrt(np.sum((D ** 2), axis=0)) + # compute Gramian (for efficiency) + DTD = np.dot(D.T,D) + + successCounter = 0 + trial = numTrials + while trial > 0: + # generate sparse signal with a minimum non-zero value + x = np.zeros((ambientDimension, 1)) + idx2 = idx + generator.shuffle(idx2) + idx3 = idx2[:sparsity] + while np.min(np.abs(x[idx3,0])) < 1e-10 : + x[idx3,0] = generator.randn(sparsity) + # sense sparse signal + y = np.dot(D, x) + + # Nic: Use sparsify OMP function (translated from Matlab) + ompopts = dict({'stopCrit':'M', 'stopTol':2*sparsity}) + starttime = time.time() # start timer + x_r2, errs, times = greed_omp_qr(y.squeeze().copy(), D.copy(), D.shape[1], ompopts) + t2all = t2all + time.time() - starttime # stop timer + idx_r2 = np.nonzero(x_r2)[0] + + # run to two times expected sparsity, or tolerance + # why? Often times, OMP can retrieve the correct solution + # when it is run for more than the expected sparsity + #x_r, idx_r = omp_qr(y,D,DTD,2*sparsity,1e-5) + # Nic: adjust tolerance to match with other function + starttime = time.time() # start timer + x_r, idx_r = omp_qr(y.copy(),D.copy(),DTD.copy(),2*sparsity,numMeasurements*1e-14/np.vdot(y,y)) + t1all = t1all + time.time() - starttime # stop timer + + # Nic: test sklearn omp + starttime = time.time() # start timer + x_r3 = omp.orthogonal_mp(D.copy(), y.copy(), 2*sparsity, tol=numMeasurements*1e-14/np.vdot(y,y), precompute_gram=False, copy_X=True) + idx_r3 = np.nonzero(x_r3)[0] + t3all = t3all + time.time() - starttime # stop timer + + # Nic: compare results + print 'diff1 = ',np.linalg.norm(x_r.squeeze() - x_r2.squeeze()) + print 'diff2 = ',np.linalg.norm(x_r.squeeze() - x_r3.squeeze()) + print 'diff3 = ',np.linalg.norm(x_r2.squeeze() - x_r3.squeeze()) + print "Bob's total time = ", t1all + print "Nic's total time = ", t2all + print "Skl's total time = ", t3all + if np.linalg.norm(x_r.squeeze() - x_r2.squeeze()) > 1e-10 or \ + np.linalg.norm(x_r.squeeze() - x_r3.squeeze()) > 1e-10 or \ + np.linalg.norm(x_r2.squeeze() - x_r3.squeeze()) > 1e-10: + print "STOP: Different results" + print "Bob's residual: ||y - D x_r ||_2 = ",np.linalg.norm(y.squeeze() - np.dot(D,x_r).squeeze()) + print "Nic's residual: ||y - D x_r ||_2 = ",np.linalg.norm(y.squeeze() - np.dot(D,x_r2).squeeze()) + print "Skl's residual: ||y - D x_r ||_2 = ",np.linalg.norm(y.squeeze() - np.dot(D,x_r3).squeeze()) + #raise ValueError("Different results") + + # debais to remove small entries + for nn in idx_r: + if abs(x_r[nn]) < 1e-10: + x_r[nn] = 0 + + # exact recovery condition using support + #if sorted(np.flatnonzero(x_r)) == sorted(np.flatnonzero(x)): + # successCounter += 1 + # exact recovery condition using error in solution + error = x - x_r + """ the following is the exact recovery condition in: A. Maleki + and D. L. Donoho, "Optimally tuned iterative reconstruction + algorithms for compressed sensing," IEEE J. Selected Topics + in Signal Process., vol. 4, pp. 330-341, Apr. 2010. """ + if np.vdot(error,error) < np.vdot(x,x)*1e-4: + successCounter += 1 + trial -= 1 + + success[rhoCounter,deltaCounter] = successCounter + if successCounter == 0: + break + + rhoCounter += 1 + #np.savetxt('test.txt',success,fmt='#2.1d',delimiter=',') + deltaCounter += 1 + + +def omp_qr(x, dict, D, natom, tolerance): + """ Recover x using QR implementation of OMP + + Parameter + --------- + x: measurements + dict: dictionary + D: Gramian of dictionary + natom: iterations + tolerance: error tolerance + + Return + ------ + x_hat : estimate of x + gamma : indices where non-zero + + For more information, see http://media.aau.dk/null_space_pursuits/2011/10/efficient-omp.html + """ + msize, dictsize = dict.shape + normr2 = np.vdot(x,x) + normtol2 = tolerance*normr2 + R = np.zeros((natom,natom)) + Q = np.zeros((msize,natom)) + gamma = [] + + # find initial projections + origprojections = np.dot(x.T,dict) + origprojectionsT = origprojections.T + projections = origprojections.copy(); + + k = 0 + while (normr2 > normtol2) and (k < natom): + # find index of maximum magnitude projection + newgam = np.argmax(np.abs(projections ** 2)) + gamma.append(newgam) + # update QR factorization, projections, and residual energy + if k == 0: + R[0,0] = 1 + Q[:,0] = dict[:,newgam].copy() + # update projections + QtempQtempT = np.outer(Q[:,0],Q[:,0]) + projections -= np.dot(x.T, np.dot(QtempQtempT,dict)) + # update residual energy + normr2 -= np.vdot(x, np.dot(QtempQtempT,x)) + else: + w = scipy.linalg.solve_triangular(R[0:k,0:k],D[gamma[0:k],newgam],trans=1) + R[k,k] = np.sqrt(1-np.vdot(w,w)) + R[0:k,k] = w.copy() + Q[:,k] = (dict[:,newgam] - np.dot(QtempQtempT,dict[:,newgam]))/R[k,k] + QkQkT = np.outer(Q[:,k],Q[:,k]) + xTQkQkT = np.dot(x.T,QkQkT) + QtempQtempT += QkQkT + # update projections + projections -= np.dot(xTQkQkT,dict) + # update residual energy + normr2 -= np.dot(xTQkQkT,x) + + k += 1 + + # build solution + tempR = R[0:k,0:k] + w = scipy.linalg.solve_triangular(tempR,origprojectionsT[gamma[0:k]],trans=1) + x_hat = np.zeros((dictsize,1)) + x_hat[gamma[0:k]] = scipy.linalg.solve_triangular(tempR,w) + + return x_hat, gamma + + +#function [s, err_mse, iter_time]=greed_omp_qr(x,A,m,varargin) +def greed_omp_qr(x,A,m,opts=[]): +# greed_omp_qr: Orthogonal Matching Pursuit algorithm based on QR +# factorisation +# Nic: translated to Python on 19.10.2011. Original Matlab Code by Thomas Blumensath +########################################################################### +# Usage +# [s, err_mse, iter_time]=greed_omp_qr(x,P,m,'option_name','option_value') +########################################################################### +########################################################################### +# Input +# Mandatory: +# x Observation vector to be decomposed +# P Either: +# 1) An nxm matrix (n must be dimension of x) +# 2) A function handle (type "help function_format" +# for more information) +# Also requires specification of P_trans option. +# 3) An object handle (type "help object_format" for +# more information) +# m length of s +# +# Possible additional options: +# (specify as many as you want using 'option_name','option_value' pairs) +# See below for explanation of options: +#__________________________________________________________________________ +# option_name | available option_values | default +#-------------------------------------------------------------------------- +# stopCrit | M, corr, mse, mse_change | M +# stopTol | number (see below) | n/4 +# P_trans | function_handle (see below) | +# maxIter | positive integer (see below) | n +# verbose | true, false | false +# start_val | vector of length m | zeros +# +# Available stopping criteria : +# M - Extracts exactly M = stopTol elements. +# corr - Stops when maximum correlation between +# residual and atoms is below stopTol value. +# mse - Stops when mean squared error of residual +# is below stopTol value. +# mse_change - Stops when the change in the mean squared +# error falls below stopTol value. +# +# stopTol: Value for stopping criterion. +# +# P_trans: If P is a function handle, then P_trans has to be specified and +# must be a function handle. +# +# maxIter: Maximum number of allowed iterations. +# +# verbose: Logical value to allow algorithm progress to be displayed. +# +# start_val: Allows algorithms to start from partial solution. +# +########################################################################### +# Outputs +# s Solution vector +# err_mse Vector containing mse of approximation error for each +# iteration +# iter_time Vector containing computation times for each iteration +# +########################################################################### +# Description +# greed_omp_qr performs a greedy signal decomposition. +# In each iteration a new element is selected depending on the inner +# product between the current residual and columns in P. +# The non-zero elements of s are approximated by orthogonally projecting +# x onto the selected elements in each iteration. +# The algorithm uses QR decomposition. +# +# See Also +# greed_omp_chol, greed_omp_cg, greed_omp_cgp, greed_omp_pinv, +# greed_omp_linsolve, greed_gp, greed_nomp +# +# Copyright (c) 2007 Thomas Blumensath +# +# The University of Edinburgh +# Email: thomas.blumensath@ed.ac.uk +# Comments and bug reports welcome +# +# This file is part of sparsity Version 0.1 +# Created: April 2007 +# +# Part of this toolbox was developed with the support of EPSRC Grant +# D000246/1 +# +# Please read COPYRIGHT.m for terms and conditions. + + ########################################################################### + # Default values and initialisation + ########################################################################### + #[n1 n2]=size(x); + #n1,n2 = x.shape + #if n2 == 1 + # n=n1; + #elseif n1 == 1 + # x=x'; + # n=n2; + #else + # display('x must be a vector.'); + # return + #end + if x.ndim != 1: + print 'x must be a vector.' + return + n = x.size + + #sigsize = x'*x/n; + sigsize = np.vdot(x,x)/n; + initial_given = 0; + err_mse = np.array([]); + iter_time = np.array([]); + STOPCRIT = 'M'; + STOPTOL = math.ceil(n/4.0); + MAXITER = n; + verbose = False; + s_initial = np.zeros(m); + + if verbose: + print 'Initialising...' + #end + + ########################################################################### + # Output variables + ########################################################################### + #switch nargout + # case 3 + # comp_err=true; + # comp_time=true; + # case 2 + # comp_err=true; + # comp_time=false; + # case 1 + # comp_err=false; + # comp_time=false; + # case 0 + # error('Please assign output variable.') + # otherwise + # error('Too many output arguments specified') + #end + if 'nargout' in opts: + if opts['nargout'] == 3: + comp_err = True + comp_time = True + elif opts['nargout'] == 2: + comp_err = True + comp_time = False + elif opts['nargout'] == 1: + comp_err = False + comp_time = False + elif opts['nargout'] == 0: + print 'Please assign output variable.' + return + else: + print 'Too many output arguments specified' + return + else: + # If not given, make default nargout = 3 + # and add nargout to options + opts['nargout'] = 3 + comp_err = True + comp_time = True + + ########################################################################### + # Look through options + ########################################################################### + # Put option into nice format + #Options={}; + #OS=nargin-3; + #c=1; + #for i=1:OS + # if isa(varargin{i},'cell') + # CellSize=length(varargin{i}); + # ThisCell=varargin{i}; + # for j=1:CellSize + # Options{c}=ThisCell{j}; + # c=c+1; + # end + # else + # Options{c}=varargin{i}; + # c=c+1; + # end + #end + #OS=length(Options); + #if rem(OS,2) + # error('Something is wrong with argument name and argument value pairs.') + #end + # + #for i=1:2:OS + # switch Options{i} + # case {'stopCrit'} + # if (strmatch(Options{i+1},{'M'; 'corr'; 'mse'; 'mse_change'},'exact')); + # STOPCRIT = Options{i+1}; + # else error('stopCrit must be char string [M, corr, mse, mse_change]. Exiting.'); end + # case {'stopTol'} + # if isa(Options{i+1},'numeric') ; STOPTOL = Options{i+1}; + # else error('stopTol must be number. Exiting.'); end + # case {'P_trans'} + # if isa(Options{i+1},'function_handle'); Pt = Options{i+1}; + # else error('P_trans must be function _handle. Exiting.'); end + # case {'maxIter'} + # if isa(Options{i+1},'numeric'); MAXITER = Options{i+1}; + # else error('maxIter must be a number. Exiting.'); end + # case {'verbose'} + # if isa(Options{i+1},'logical'); verbose = Options{i+1}; + # else error('verbose must be a logical. Exiting.'); end + # case {'start_val'} + # if isa(Options{i+1},'numeric') & length(Options{i+1}) == m ; + # s_initial = Options{i+1}; + # initial_given=1; + # else error('start_val must be a vector of length m. Exiting.'); end + # otherwise + # error('Unrecognised option. Exiting.') + # end + #end + if 'stopCrit' in opts: + STOPCRIT = opts['stopCrit'] + if 'stopTol' in opts: + if hasattr(opts['stopTol'], '__int__'): # check if numeric + STOPTOL = opts['stopTol'] + else: + raise TypeError('stopTol must be number. Exiting.') + if 'P_trans' in opts: + if hasattr(opts['P_trans'], '__call__'): # check if function handle + Pt = opts['P_trans'] + else: + raise TypeError('P_trans must be function _handle. Exiting.') + if 'maxIter' in opts: + if hasattr(opts['maxIter'], '__int__'): # check if numeric + MAXITER = opts['maxIter'] + else: + raise TypeError('maxIter must be a number. Exiting.') + if 'verbose' in opts: + # TODO: Should check here if is logical + verbose = opts['verbose'] + if 'start_val' in opts: + # TODO: Should check here if is numeric + if opts['start_val'].size == m: + s_initial = opts['start_val'] + initial_given = 1 + else: + raise ValueError('start_val must be a vector of length m. Exiting.') + # Don't exit if unknown option is given, simply ignore it + + #if strcmp(STOPCRIT,'M') + # maxM=STOPTOL; + #else + # maxM=MAXITER; + #end + if STOPCRIT == 'M': + maxM = STOPTOL + else: + maxM = MAXITER + + # if nargout >=2 + # err_mse = zeros(maxM,1); + # end + # if nargout ==3 + # iter_time = zeros(maxM,1); + # end + if opts['nargout'] >= 2: + err_mse = np.zeros(maxM) + if opts['nargout'] == 3: + iter_time = np.zeros(maxM) + + ########################################################################### + # Make P and Pt functions + ########################################################################### + #if isa(A,'float') P =@(z) A*z; Pt =@(z) A'*z; + #elseif isobject(A) P =@(z) A*z; Pt =@(z) A'*z; + #elseif isa(A,'function_handle') + # try + # if isa(Pt,'function_handle'); P=A; + # else error('If P is a function handle, Pt also needs to be a function handle. Exiting.'); end + # catch error('If P is a function handle, Pt needs to be specified. Exiting.'); end + #else error('P is of unsupported type. Use matrix, function_handle or object. Exiting.'); end + if hasattr(A, '__call__'): + if hasattr(Pt, '__call__'): + P = A + else: + raise TypeError('If P is a function handle, Pt also needs to be a function handle.') + else: + # TODO: should check here if A is matrix + P = lambda z: np.dot(A,z) + Pt = lambda z: np.dot(A.T,z) + + ########################################################################### + # Random Check to see if dictionary is normalised + ########################################################################### + # mask=zeros(m,1); + # mask(ceil(rand*m))=1; + # nP=norm(P(mask)); + # if abs(1-nP)>1e-3; + # display('Dictionary appears not to have unit norm columns.') + # end + mask = np.zeros(m) + mask[math.floor(np.random.rand() * m)] = 1 + nP = np.linalg.norm(P(mask)) + if abs(1-nP) > 1e-3: + print 'Dictionary appears not to have unit norm columns.' + #end + + ########################################################################### + # Check if we have enough memory and initialise + ########################################################################### + # try Q=zeros(n,maxM); + # catch error('Variable size is too large. Please try greed_omp_chol algorithm or reduce MAXITER.') + # end + # try R=zeros(maxM); + # catch error('Variable size is too large. Please try greed_omp_chol algorithm or reduce MAXITER.') + # end + try: + Q = np.zeros((n,maxM)) + except: + print 'Variable size is too large. Please try greed_omp_chol algorithm or reduce MAXITER.' + raise + try: + R = np.zeros((maxM, maxM)) + except: + print 'Variable size is too large. Please try greed_omp_chol algorithm or reduce MAXITER.' + raise + + ########################################################################### + # Do we start from zero or not? + ########################################################################### + #if initial_given ==1; + # IN = find(s_initial); + # if ~isempty(IN) + # Residual = x-P(s_initial); + # lengthIN=length(IN); + # z=[]; + # for k=1:length(IN) + # # Extract new element + # mask=zeros(m,1); + # mask(IN(k))=1; + # new_element=P(mask); + # + # # Orthogonalise new element + # qP=Q(:,1:k-1)'*new_element; + # q=new_element-Q(:,1:k-1)*(qP); + # + # nq=norm(q); + # q=q/nq; + # # Update QR factorisation + # R(1:k-1,k)=qP; + # R(k,k)=nq; + # Q(:,k)=q; + # + # z(k)=q'*x; + # end + # s = s_initial; + # Residual=x-Q(:,k)*z; + # oldERR = Residual'*Residual/n; + # else + # IN = []; + # Residual = x; + # s = s_initial; + # sigsize = x'*x/n; + # oldERR = sigsize; + # k=0; + # z=[]; + # end + # + #else + # IN = []; + # Residual = x; + # s = s_initial; + # sigsize = x'*x/n; + # oldERR = sigsize; + # k=0; + # z=[]; + #end + if initial_given == 1: + #IN = find(s_initial); + IN = np.nonzero(s_initial)[0].tolist() + #if ~isempty(IN) + if IN.size > 0: + Residual = x - P(s_initial) + lengthIN = IN.size + z = np.array([]) + #for k=1:length(IN) + for k in np.arange(IN.size): + # Extract new element + mask = np.zeros(m) + mask[IN[k]] = 1 + new_element = P(mask) + + # Orthogonalise new element + #qP=Q(:,1:k-1)'*new_element; + if k-1 >= 0: + qP = np.dot(Q[:,0:k].T , new_element) + #q=new_element-Q(:,1:k-1)*(qP); + q = new_element - np.dot(Q[:,0:k] , qP) + + nq = np.linalg.norm(q) + q = q / nq + # Update QR factorisation + R[0:k,k] = qP + R[k,k] = nq + Q[:,k] = q + else: + q = new_element + + nq = np.linalg.norm(q) + q = q / nq + # Update QR factorisation + R[k,k] = nq + Q[:,k] = q + + z[k] = np.dot(q.T , x) + #end + s = s_initial.copy() + Residual = x - np.dot(Q[:,k] , z) + oldERR = np.vdot(Residual , Residual) / n; + else: + #IN = np.array([], dtype = int) + IN = np.array([], dtype = int).tolist() + Residual = x.copy() + s = s_initial.copy() + sigsize = np.vdot(x , x) / n + oldERR = sigsize + k = 0 + #z = np.array([]) + z = [] + #end + + else: + #IN = np.array([], dtype = int) + IN = np.array([], dtype = int).tolist() + Residual = x.copy() + s = s_initial.copy() + sigsize = np.vdot(x , x) / n + oldERR = sigsize + k = 0 + #z = np.array([]) + z = [] + #end + + ########################################################################### + # Main algorithm + ########################################################################### + # if verbose + # display('Main iterations...') + # end + # tic + # t=0; + # DR=Pt(Residual); + # done = 0; + # iter=1; + if verbose: + print 'Main iterations...' + tic = time.time() + t = 0 + DR = Pt(Residual) + done = 0 + iter = 1 + + #while ~done + # + # # Select new element + # DR(IN)=0; + # # Nic: replace selection with random variable + # # i.e. Randomized OMP!! + # # DON'T FORGET ABOUT THIS!! + # [v I]=max(abs(DR)); + # #I = randp(exp(abs(DR).^2 ./ (norms.^2)'), [1 1]); + # IN=[IN I]; + # + # + # k=k+1; + # # Extract new element + # mask=zeros(m,1); + # mask(IN(k))=1; + # new_element=P(mask); + # + # # Orthogonalise new element + # qP=Q(:,1:k-1)'*new_element; + # q=new_element-Q(:,1:k-1)*(qP); + # + # nq=norm(q); + # q=q/nq; + # # Update QR factorisation + # R(1:k-1,k)=qP; + # R(k,k)=nq; + # Q(:,k)=q; + # + # z(k)=q'*x; + # + # # New residual + # Residual=Residual-q*(z(k)); + # DR=Pt(Residual); + # + # ERR=Residual'*Residual/n; + # if comp_err + # err_mse(iter)=ERR; + # end + # + # if comp_time + # iter_time(iter)=toc; + # end + # + ############################################################################ + ## Are we done yet? + ############################################################################ + # + # if strcmp(STOPCRIT,'M') + # if iter >= STOPTOL + # done =1; + # elseif verbose && toc-t>10 + # display(sprintf('Iteration #i. --- #i iterations to go',iter ,STOPTOL-iter)) + # t=toc; + # end + # elseif strcmp(STOPCRIT,'mse') + # if comp_err + # if err_mse(iter)<STOPTOL; + # done = 1; + # elseif verbose && toc-t>10 + # display(sprintf('Iteration #i. --- #i mse',iter ,err_mse(iter))) + # t=toc; + # end + # else + # if ERR<STOPTOL; + # done = 1; + # elseif verbose && toc-t>10 + # display(sprintf('Iteration #i. --- #i mse',iter ,ERR)) + # t=toc; + # end + # end + # elseif strcmp(STOPCRIT,'mse_change') && iter >=2 + # if comp_err && iter >=2 + # if ((err_mse(iter-1)-err_mse(iter))/sigsize <STOPTOL); + # done = 1; + # elseif verbose && toc-t>10 + # display(sprintf('Iteration #i. --- #i mse change',iter ,(err_mse(iter-1)-err_mse(iter))/sigsize )) + # t=toc; + # end + # else + # if ((oldERR - ERR)/sigsize < STOPTOL); + # done = 1; + # elseif verbose && toc-t>10 + # display(sprintf('Iteration #i. --- #i mse change',iter ,(oldERR - ERR)/sigsize)) + # t=toc; + # end + # end + # elseif strcmp(STOPCRIT,'corr') + # if max(abs(DR)) < STOPTOL; + # done = 1; + # elseif verbose && toc-t>10 + # display(sprintf('Iteration #i. --- #i corr',iter ,max(abs(DR)))) + # t=toc; + # end + # end + # + # # Also stop if residual gets too small or maxIter reached + # if comp_err + # if err_mse(iter)<1e-16 + # display('Stopping. Exact signal representation found!') + # done=1; + # end + # else + # + # + # if iter>1 + # if ERR<1e-16 + # display('Stopping. Exact signal representation found!') + # done=1; + # end + # end + # end + # + # if iter >= MAXITER + # display('Stopping. Maximum number of iterations reached!') + # done = 1; + # end + # + ############################################################################ + ## If not done, take another round + ############################################################################ + # + # if ~done + # iter=iter+1; + # oldERR=ERR; + # end + #end + while not done: + + # Select new element + DR[IN]=0 + #[v I]=max(abs(DR)); + #v = np.abs(DR).max() + I = np.abs(DR).argmax() + #IN = np.concatenate((IN,I)) + IN.append(I) + + + #k = k + 1 Move to end, since is zero based + + # Extract new element + mask = np.zeros(m) + mask[IN[k]] = 1 + new_element = P(mask) + + # Orthogonalise new element + if k-1 >= 0: + qP = np.dot(Q[:,0:k].T , new_element) + q = new_element - np.dot(Q[:,0:k] , qP) + + nq = np.linalg.norm(q) + q = q/nq + # Update QR factorisation + R[0:k,k] = qP + R[k,k] = nq + Q[:,k] = q + else: + q = new_element + + nq = np.linalg.norm(q) + q = q/nq + # Update QR factorisation + R[k,k] = nq + Q[:,k] = q + + #z[k]=np.vdot(q , x) + z.append(np.vdot(q , x)) + + # New residual + Residual = Residual - q * (z[k]) + DR = Pt(Residual) + + ERR = np.vdot(Residual , Residual) / n + if comp_err: + err_mse[iter-1] = ERR + #end + + if comp_time: + iter_time[iter-1] = time.time() - tic + #end + + ########################################################################### + # Are we done yet? + ########################################################################### + if STOPCRIT == 'M': + if iter >= STOPTOL: + done = 1 + elif verbose and time.time() - t > 10.0/1000: # time() returns sec + #display(sprintf('Iteration #i. --- #i iterations to go',iter ,STOPTOL-iter)) + print 'Iteration '+iter+'. --- '+(STOPTOL-iter)+' iterations to go' + t = time.time() + #end + elif STOPCRIT =='mse': + if comp_err: + if err_mse[iter-1] < STOPTOL: + done = 1 + elif verbose and time.time() - t > 10.0/1000: # time() returns sec + #display(sprintf('Iteration #i. --- #i mse',iter ,err_mse(iter))) + print 'Iteration '+iter+'. --- '+err_mse[iter-1]+' mse' + t = time.time() + #end + else: + if ERR < STOPTOL: + done = 1 + elif verbose and time.time() - t > 10.0/1000: # time() returns sec + #display(sprintf('Iteration #i. --- #i mse',iter ,ERR)) + print 'Iteration '+iter+'. --- '+ERR+' mse' + t = time.time() + #end + #end + elif STOPCRIT == 'mse_change' and iter >=2: + if comp_err and iter >=2: + if ((err_mse[iter-2] - err_mse[iter-1])/sigsize < STOPTOL): + done = 1 + elif verbose and time.time() - t > 10.0/1000: # time() returns sec + #display(sprintf('Iteration #i. --- #i mse change',iter ,(err_mse(iter-1)-err_mse(iter))/sigsize )) + print 'Iteration '+iter+'. --- '+((err_mse[iter-2]-err_mse[iter-1])/sigsize)+' mse change' + t = time.time() + #end + else: + if ((oldERR - ERR)/sigsize < STOPTOL): + done = 1 + elif verbose and time.time() - t > 10.0/1000: # time() returns sec + #display(sprintf('Iteration #i. --- #i mse change',iter ,(oldERR - ERR)/sigsize)) + print 'Iteration '+iter+'. --- '+((oldERR - ERR)/sigsize)+' mse change' + t = time.time() + #end + #end + elif STOPCRIT == 'corr': + if np.abs(DR).max() < STOPTOL: + done = 1 + elif verbose and time.time() - t > 10.0/1000: # time() returns sec + #display(sprintf('Iteration #i. --- #i corr',iter ,max(abs(DR)))) + print 'Iteration '+iter+'. --- '+(np.abs(DR).max())+' corr' + t = time.time() + #end + #end + + # Also stop if residual gets too small or maxIter reached + if comp_err: + if err_mse[iter-1] < 1e-14: + done = 1 + # Nic: added verbose check + if verbose: + print 'Stopping. Exact signal representation found!' + #end + else: + if iter > 1: + if ERR < 1e-14: + done = 1 + # Nic: added verbose check + if verbose: + print 'Stopping. Exact signal representation found!' + #end + #end + #end + + + if iter >= MAXITER: + done = 1 + # Nic: added verbose check + if verbose: + print 'Stopping. Maximum number of iterations reached!' + #end + + ########################################################################### + # If not done, take another round + ########################################################################### + if not done: + iter = iter + 1 + oldERR = ERR + #end + + # Moved here from front, since we are 0-based + k = k + 1 + #end + + ########################################################################### + # Now we can solve for s by back-substitution + ########################################################################### + #s(IN)=R(1:k,1:k)\z(1:k)'; + s[IN] = scipy.linalg.solve(R[0:k,0:k] , np.array(z[0:k])) + + ########################################################################### + # Only return as many elements as iterations + ########################################################################### + if opts['nargout'] >= 2: + err_mse = err_mse[0:iter-1] + #end + if opts['nargout'] == 3: + iter_time = iter_time[0:iter-1] + #end + if verbose: + print 'Done' + #end + + # Return + if opts['nargout'] == 1: + return s + elif opts['nargout'] == 2: + return s, err_mse + elif opts['nargout'] == 3: + return s, err_mse, iter_time + + # Change history + # + # 8 of Februray: Algo does no longer stop if dictionary is not normaliesd. + + # End of greed_omp_qr() function + #-------------------------------- + +if __name__ == '__main__': + print ('Running problem suite') + ambientDimension = 400 + numGradations = 30 + numTrials = 1 + + #import cProfile + #cProfile.run('runProblemSuite(ambientDimension,numGradations,numTrials)','profres') + runProblemSuite(ambientDimension,numGradations,numTrials) + print "Done" + + #import pstats + #p = pstats.Stats('D:\Nic\Dev2\profres') + #p.sort_stats('cumulative').print_stats(10)