diff options
Diffstat (limited to 'Optimization')
| -rw-r--r-- | Optimization/cgsolve.m | 76 | ||||
| -rw-r--r-- | Optimization/l1dantzig_pd.m | 234 | ||||
| -rw-r--r-- | Optimization/l1decode_pd.m | 175 | ||||
| -rw-r--r-- | Optimization/l1eq_pd.m | 209 | ||||
| -rw-r--r-- | Optimization/l1qc_logbarrier.m | 116 | ||||
| -rw-r--r-- | Optimization/l1qc_newton.m | 150 | ||||
| -rw-r--r-- | Optimization/tvdantzig_logbarrier.m | 120 | ||||
| -rw-r--r-- | Optimization/tvdantzig_newton.m | 185 | ||||
| -rw-r--r-- | Optimization/tveq_logbarrier.m | 118 | ||||
| -rw-r--r-- | Optimization/tveq_newton.m | 180 | ||||
| -rw-r--r-- | Optimization/tvqc_logbarrier.m | 121 | ||||
| -rw-r--r-- | Optimization/tvqc_newton.m | 176 |
12 files changed, 1860 insertions, 0 deletions
diff --git a/Optimization/cgsolve.m b/Optimization/cgsolve.m new file mode 100644 index 0000000..dbbc4d8 --- /dev/null +++ b/Optimization/cgsolve.m @@ -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; + diff --git a/Optimization/l1dantzig_pd.m b/Optimization/l1dantzig_pd.m new file mode 100644 index 0000000..6a57dea --- /dev/null +++ b/Optimization/l1dantzig_pd.m @@ -0,0 +1,234 @@ +% l1dantzig_pd.m +% +% Solves +% min_x ||x||_1 subject to ||A'(Ax-b)||_\infty <= epsilon +% +% Recast as linear program +% min_{x,u} sum(u) s.t. x - u <= 0 +% -x - u <= 0 +% A'(Ax-b) - epsilon <= 0 +% -A'(Ax-b) - epsilon <= 0 +% and use primal-dual interior point method. +% +% Usage: xp = l1dantzig_pd(x0, A, At, b, epsilon, pdtol, pdmaxiter, 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 or Nx1 vector of correlation constraints +% +% pdtol - Tolerance for primal-dual algorithm (algorithm terminates if +% the duality gap is less than pdtol). +% Default = 1e-3. +% +% pdmaxiter - Maximum number of primal-dual 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 = l1dantzig_pd(x0, A, At, b, epsilon, pdtol, pdmaxiter, cgtol, cgmaxiter) + +largescale = isa(A,'function_handle'); + +if (nargin < 6), pdtol = 1e-3; end +if (nargin < 7), pdmaxiter = 50; end +if (nargin < 8), cgtol = 1e-8; end +if (nargin < 9), cgmaxiter = 200; end + +N = length(x0); + +alpha = 0.01; +beta = 0.5; +mu = 10; + +gradf0 = [zeros(N,1); ones(N,1)]; + + +% starting point --- make sure that it is feasible +if (largescale) + if (max( abs(At(A(x0) - b)) - epsilon ) > 0) + disp('Starting point infeasible; using x0 = At*inv(AAt)*y.'); + AAt = @(z) A(At(z)); + [w, cgres] = 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 (max(abs(A'*(A*x0 - b)) - epsilon ) > 0) + 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)); + +% set up for the first iteration +if (largescale) + Atr = At(A(x) - b); +else + Atr = A'*(A*x - b); +end +fu1 = x - u; +fu2 = -x - u; +fe1 = Atr - epsilon; +fe2 = -Atr - epsilon; +lamu1 = -(1./fu1); +lamu2 = -(1./fu2); +lame1 = -(1./fe1); +lame2 = -(1./fe2); +if (largescale) + AtAv = At(A(lame1-lame2)); +else + AtAv = A'*(A*(lame1-lame2)); +end + +% sdg = surrogate duality gap +sdg = -[fu1; fu2; fe1; fe2]'*[lamu1; lamu2; lame1; lame2]; +tau = mu*(4*N)/sdg; + +% residuals +rdual = gradf0 + [lamu1-lamu2 + AtAv; -lamu1-lamu2]; +rcent = -[lamu1.*fu1; lamu2.*fu2; lame1.*fe1; lame2.*fe2] - (1/tau); +resnorm = norm([rdual; rcent]); + +% iterations +pditer = 0; +done = (sdg < pdtol) | (pditer >= pdmaxiter); +while (~done) + + % solve for step direction + w2 = - 1 - (1/tau)*(1./fu1 + 1./fu2); + + sig11 = -lamu1./fu1 - lamu2./fu2; + sig12 = lamu1./fu1 - lamu2./fu2; + siga = -(lame1./fe1 + lame2./fe2); + sigx = sig11 - sig12.^2./sig11; + + if (largescale) + w1 = -(1/tau)*( At(A(1./fe2-1./fe1)) + 1./fu2 - 1./fu1 ); + w1p = w1 - (sig12./sig11).*w2; + hpfun = @(z) At(A(siga.*At(A(z)))) + sigx.*z; + [dx, cgres, cgiter] = cgsolve(hpfun, 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; + return + end + AtAdx = At(A(dx)); + else + w1 = -(1/tau)*( A'*(A*(1./fe2-1./fe1)) + 1./fu2 - 1./fu1 ); + w1p = w1 - (sig12./sig11).*w2; + Hp = A'*(A*sparse(diag(siga))*A')*A + diag(sigx); + opts.POSDEF = true; opts.SYM = true; + [dx, hcond] = linsolve(Hp, w1p,opts); + if (hcond < 1e-14) + disp('Matrix ill-conditioned. Returning previous iterate. (See Section 4 of notes for more information.)'); + xp = x; + return + end + AtAdx = A'*(A*dx); + end + du = w2./sig11 - (sig12./sig11).*dx; + + dlamu1 = -(lamu1./fu1).*(dx-du) - lamu1 - (1/tau)*1./fu1; + dlamu2 = -(lamu2./fu2).*(-dx-du) - lamu2 - (1/tau)*1./fu2; + dlame1 = -(lame1./fe1).*(AtAdx) - lame1 - (1/tau)*1./fe1; + dlame2 = -(lame2./fe2).*(-AtAdx) - lame2 - (1/tau)*1./fe2; + if (largescale) + AtAdv = At(A(dlame1-dlame2)); + else + AtAdv = A'*(A*(dlame1-dlame2)); + end + + + % find minimal step size that keeps ineq functions < 0, dual vars > 0 + iu1 = find(dlamu1 < 0); iu2 = find(dlamu2 < 0); + ie1 = find(dlame1 < 0); ie2 = find(dlame2 < 0); + ifu1 = find((dx-du) > 0); ifu2 = find((-dx-du) > 0); + ife1 = find(AtAdx > 0); ife2 = find(-AtAdx > 0); + smax = min(1,min([... + -lamu1(iu1)./dlamu1(iu1); -lamu2(iu2)./dlamu2(iu2); ... + -lame1(ie1)./dlame1(ie1); -lame2(ie2)./dlame2(ie2); ... + -fu1(ifu1)./(dx(ifu1)-du(ifu1)); -fu2(ifu2)./(-dx(ifu2)-du(ifu2)); ... + -fe1(ife1)./AtAdx(ife1); -fe2(ife2)./(-AtAdx(ife2)) ])); + s = 0.99*smax; + + % backtracking line search + suffdec = 0; + backiter = 0; + while (~suffdec) + xp = x + s*dx; up = u + s*du; + Atrp = Atr + s*AtAdx; AtAvp = AtAv + s*AtAdv; + fu1p = fu1 + s*(dx-du); fu2p = fu2 + s*(-dx-du); + fe1p = fe1 + s*AtAdx; fe2p = fe2 + s*(-AtAdx); + lamu1p = lamu1 + s*dlamu1; lamu2p = lamu2 + s*dlamu2; + lame1p = lame1 + s*dlame1; lame2p = lame2 + s*dlame2; + rdp = gradf0 + [lamu1p-lamu2p + AtAvp; -lamu1p-lamu2p]; + rcp = -[lamu1p.*fu1p; lamu2p.*fu2p; lame1p.*fe1p; lame2p.*fe2p] - (1/tau); + suffdec = (norm([rdp; rcp]) <= (1-alpha*s)*resnorm); + s = beta*s; + backiter = backiter+1; + if (backiter > 32) + disp('Stuck backtracking, returning last iterate. (See Section 4 of notes for more information.)') + xp = x; + return + end + end + + % setup for next iteration + x = xp; u = up; + Atr = Atrp; AtAv = AtAvp; + fu1 = fu1p; fu2 = fu2p; + fe1 = fe1p; fe2 = fe2p; + lamu1 = lamu1p; lamu2 = lamu2p; + lame1 = lame1p; lame2 = lame2p; + + sdg = -[fu1; fu2; fe1; fe2]'*[lamu1; lamu2; lame1; lame2]; + tau = mu*(4*N)/sdg; + + rdual = rdp; + rcent = -[lamu1.*fu1; lamu2.*fu2; lame1.*fe1; lame2.*fe2] - (1/tau); + resnorm = norm([rdual; rcent]); + + pditer = pditer+1; + done = (sdg < pdtol) | (pditer >= pdmaxiter); + + disp(sprintf('Iteration = %d, tau = %8.3e, Primal = %8.3e, PDGap = %8.3e, Dual res = %8.3e',... + pditer, tau, sum(u), sdg, norm(rdual))); + if (largescale) + disp(sprintf(' CG Res = %8.3e, CG Iter = %d', cgres, cgiter)); + else + disp(sprintf(' H11p condition number = %8.3e', hcond)); + end + +end diff --git a/Optimization/l1decode_pd.m b/Optimization/l1decode_pd.m new file mode 100644 index 0000000..2757404 --- /dev/null +++ b/Optimization/l1decode_pd.m @@ -0,0 +1,175 @@ +% l1decode_pd.m +% +% Decoding via linear programming. +% Solve +% min_x ||b-Ax||_1 . +% +% Recast as the linear program +% min_{x,u} sum(u) s.t. -Ax - u + y <= 0 +% Ax - u - y <= 0 +% and solve using primal-dual interior point method. +% +% Usage: xp = l1decode_pd(x0, A, At, y, pdtol, pdmaxiter, cgtol, cgmaxiter) +% +% x0 - Nx1 vector, initial point. +% +% A - Either a handle to a function that takes a N vector and returns a M +% vector, or a MxN 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 an M vector and returns an N vector. +% If A is a matrix, At is ignored. +% +% y - Mx1 observed code (M > N). +% +% pdtol - Tolerance for primal-dual algorithm (algorithm terminates if +% the duality gap is less than pdtol). +% Default = 1e-3. +% +% pdmaxiter - Maximum number of primal-dual 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 = l1decode_pd(x0, A, At, y, pdtol, pdmaxiter, cgtol, cgmaxiter) + +largescale = isa(A,'function_handle'); + +if (nargin < 5), pdtol = 1e-3; end +if (nargin < 6), pdmaxiter = 50; end +if (nargin < 7), cgtol = 1e-8; end +if (nargin < 8), cgmaxiter = 200; end + +N = length(x0); +M = length(y); + +alpha = 0.01; +beta = 0.5; +mu = 10; + +gradf0 = [zeros(N,1); ones(M,1)]; + +x = x0; +if (largescale), Ax = A(x); else Ax = A*x; end +u = (0.95)*abs(y-Ax) + (0.10)*max(abs(y-Ax)); + +fu1 = Ax - y - u; +fu2 = -Ax + y - u; + +lamu1 = -1./fu1; +lamu2 = -1./fu2; + +if (largescale), Atv = At(lamu1-lamu2); else Atv = A'*(lamu1-lamu2); end + +sdg = -(fu1'*lamu1 + fu2'*lamu2); +tau = mu*2*M/sdg; + +rcent = [-lamu1.*fu1; -lamu2.*fu2] - (1/tau); +rdual = gradf0 + [Atv; -lamu1-lamu2]; +resnorm = norm([rdual; rcent]); + +pditer = 0; +done = (sdg < pdtol)| (pditer >= pdmaxiter); +while (~done) + + pditer = pditer + 1; + + w2 = -1 - 1/tau*(1./fu1 + 1./fu2); + + sig1 = -lamu1./fu1 - lamu2./fu2; + sig2 = lamu1./fu1 - lamu2./fu2; + sigx = sig1 - sig2.^2./sig1; + + if (largescale) + w1 = -1/tau*(At(-1./fu1 + 1./fu2)); + w1p = w1 - At((sig2./sig1).*w2); + h11pfun = @(z) At(sigx.*A(z)); + [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; + return + end + Adx = A(dx); + else + w1 = -1/tau*(A'*(-1./fu1 + 1./fu2)); + w1p = w1 - A'*((sig2./sig1).*w2); + H11p = A'*(sparse(diag(sigx))*A); + 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; + return + end + Adx = A*dx; + end + + du = (w2 - sig2.*Adx)./sig1; + + dlamu1 = -(lamu1./fu1).*(Adx-du) - lamu1 - (1/tau)*1./fu1; + dlamu2 = (lamu2./fu2).*(Adx + du) -lamu2 - (1/tau)*1./fu2; + if (largescale), Atdv = At(dlamu1-dlamu2); else Atdv = A'*(dlamu1-dlamu2); end + + % make sure that the step is feasible: keeps lamu1,lamu2 > 0, fu1,fu2 < 0 + indl = find(dlamu1 < 0); indu = find(dlamu2 < 0); + s = min([1; -lamu1(indl)./dlamu1(indl); -lamu2(indu)./dlamu2(indu)]); + indl = find((Adx-du) > 0); indu = find((-Adx-du) > 0); + s = (0.99)*min([s; -fu1(indl)./(Adx(indl)-du(indl)); -fu2(indu)./(-Adx(indu)-du(indu))]); + + % backtrack + suffdec = 0; + backiter = 0; + while(~suffdec) + xp = x + s*dx; up = u + s*du; + Axp = Ax + s*Adx; Atvp = Atv + s*Atdv; + lamu1p = lamu1 + s*dlamu1; lamu2p = lamu2 + s*dlamu2; + fu1p = Axp - y - up; fu2p = -Axp + y - up; + rdp = gradf0 + [Atvp; -lamu1p-lamu2p]; + rcp = [-lamu1p.*fu1p; -lamu2p.*fu2p] - (1/tau); + suffdec = (norm([rdp; rcp]) <= (1-alpha*s)*resnorm); + s = beta*s; + backiter = backiter + 1; + if (backiter > 32) + disp('Stuck backtracking, returning last iterate. (See Section 4 of notes for more information.)') + xp = x; + return + end + end + + % next iteration + x = xp; u = up; + Ax = Axp; Atv = Atvp; + lamu1 = lamu1p; lamu2 = lamu2p; + fu1 = fu1p; fu2 = fu2p; + + % surrogate duality gap + sdg = -(fu1'*lamu1 + fu2'*lamu2); + tau = mu*2*M/sdg; + rcent = [-lamu1.*fu1; -lamu2.*fu2] - (1/tau); + rdual = rdp; + resnorm = norm([rdual; rcent]); + + done = (sdg < pdtol) | (pditer >= pdmaxiter); + + disp(sprintf('Iteration = %d, tau = %8.3e, Primal = %8.3e, PDGap = %8.3e, Dual res = %8.3e',... + pditer, tau, sum(u), sdg, norm(rdual))); + if (largescale) + disp(sprintf(' CG Res = %8.3e, CG Iter = %d', cgres, cgiter)); + else + disp(sprintf(' H11p condition number = %8.3e', hcond)); + end + +end + diff --git a/Optimization/l1eq_pd.m b/Optimization/l1eq_pd.m new file mode 100644 index 0000000..80c058e --- /dev/null +++ b/Optimization/l1eq_pd.m @@ -0,0 +1,209 @@ +% l1eq_pd.m +% +% Solve +% min_x ||x||_1 s.t. Ax = b +% +% Recast as linear program +% min_{x,u} sum(u) s.t. -u <= x <= u, Ax=b +% and use primal-dual interior point method +% +% Usage: xp = l1eq_pd(x0, A, At, b, pdtol, pdmaxiter, 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. +% +% pdtol - Tolerance for primal-dual algorithm (algorithm terminates if +% the duality gap is less than pdtol). +% Default = 1e-3. +% +% pdmaxiter - Maximum number of primal-dual 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 = l1eq_pd(x0, A, At, b, pdtol, pdmaxiter, cgtol, cgmaxiter) + +largescale = isa(A,'function_handle'); + +if (nargin < 5), pdtol = 1e-3; end +if (nargin < 6), pdmaxiter = 50; end +if (nargin < 7), cgtol = 1e-8; end +if (nargin < 8), cgmaxiter = 200; end + +N = length(x0); + +alpha = 0.01; +beta = 0.5; +mu = 10; + +gradf0 = [zeros(N,1); ones(N,1)]; + +% starting point --- make sure that it is feasible +if (largescale) + if (norm(A(x0)-b)/norm(b) > cgtol) + disp('Starting point infeasible; using x0 = At*inv(AAt)*y.'); + AAt = @(z) A(At(z)); + [w, cgres, cgiter] = 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)/norm(b) > cgtol) + 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)); + +% set up for the first iteration +fu1 = x - u; +fu2 = -x - u; +lamu1 = -1./fu1; +lamu2 = -1./fu2; +if (largescale) + v = -A(lamu1-lamu2); + Atv = At(v); + rpri = A(x) - b; +else + v = -A*(lamu1-lamu2); + Atv = A'*v; + rpri = A*x - b; +end + +sdg = -(fu1'*lamu1 + fu2'*lamu2); +tau = mu*2*N/sdg; + +rcent = [-lamu1.*fu1; -lamu2.*fu2] - (1/tau); +rdual = gradf0 + [lamu1-lamu2; -lamu1-lamu2] + [Atv; zeros(N,1)]; +resnorm = norm([rdual; rcent; rpri]); + +pditer = 0; +done = (sdg < pdtol) | (pditer >= pdmaxiter); +while (~done) + + pditer = pditer + 1; + + w1 = -1/tau*(-1./fu1 + 1./fu2) - Atv; + w2 = -1 - 1/tau*(1./fu1 + 1./fu2); + w3 = -rpri; + + sig1 = -lamu1./fu1 - lamu2./fu2; + sig2 = lamu1./fu1 - lamu2./fu2; + sigx = sig1 - sig2.^2./sig1; + + if (largescale) + w1p = w3 - A(w1./sigx - w2.*sig2./(sigx.*sig1)); + h11pfun = @(z) -A(1./sigx.*At(z)); + [dv, 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; + return + end + dx = (w1 - w2.*sig2./sig1 - At(dv))./sigx; + Adx = A(dx); + Atdv = At(dv); + else + w1p = -(w3 - A*(w1./sigx - w2.*sig2./(sigx.*sig1))); + H11p = A*(sparse(diag(1./sigx))*A'); + opts.POSDEF = true; opts.SYM = true; + [dv,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; + return + end + dx = (w1 - w2.*sig2./sig1 - A'*dv)./sigx; + Adx = A*dx; + Atdv = A'*dv; + end + + du = (w2 - sig2.*dx)./sig1; + + dlamu1 = (lamu1./fu1).*(-dx+du) - lamu1 - (1/tau)*1./fu1; + dlamu2 = (lamu2./fu2).*(dx+du) - lamu2 - 1/tau*1./fu2; + + % make sure that the step is feasible: keeps lamu1,lamu2 > 0, fu1,fu2 < 0 + indp = find(dlamu1 < 0); indn = find(dlamu2 < 0); + s = min([1; -lamu1(indp)./dlamu1(indp); -lamu2(indn)./dlamu2(indn)]); + indp = find((dx-du) > 0); indn = find((-dx-du) > 0); + s = (0.99)*min([s; -fu1(indp)./(dx(indp)-du(indp)); -fu2(indn)./(-dx(indn)-du(indn))]); + + % backtracking line search + suffdec = 0; + backiter = 0; + while (~suffdec) + xp = x + s*dx; up = u + s*du; + vp = v + s*dv; Atvp = Atv + s*Atdv; + lamu1p = lamu1 + s*dlamu1; lamu2p = lamu2 + s*dlamu2; + fu1p = xp - up; fu2p = -xp - up; + rdp = gradf0 + [lamu1p-lamu2p; -lamu1p-lamu2p] + [Atvp; zeros(N,1)]; + rcp = [-lamu1p.*fu1p; -lamu2p.*fu2p] - (1/tau); + rpp = rpri + s*Adx; + suffdec = (norm([rdp; rcp; rpp]) <= (1-alpha*s)*resnorm); + s = beta*s; + backiter = backiter + 1; + if (backiter > 32) + disp('Stuck backtracking, returning last iterate. (See Section 4 of notes for more information.)') + xp = x; + return + end + end + + + % next iteration + x = xp; u = up; + v = vp; Atv = Atvp; + lamu1 = lamu1p; lamu2 = lamu2p; + fu1 = fu1p; fu2 = fu2p; + + % surrogate duality gap + sdg = -(fu1'*lamu1 + fu2'*lamu2); + tau = mu*2*N/sdg; + rpri = rpp; + rcent = [-lamu1.*fu1; -lamu2.*fu2] - (1/tau); + rdual = gradf0 + [lamu1-lamu2; -lamu1-lamu2] + [Atv; zeros(N,1)]; + resnorm = norm([rdual; rcent; rpri]); + + done = (sdg < pdtol) | (pditer >= pdmaxiter); + + disp(sprintf('Iteration = %d, tau = %8.3e, Primal = %8.3e, PDGap = %8.3e, Dual res = %8.3e, Primal res = %8.3e',... + pditer, tau, sum(u), sdg, norm(rdual), norm(rpri))); + if (largescale) + disp(sprintf(' CG Res = %8.3e, CG Iter = %d', cgres, cgiter)); + else + disp(sprintf(' H11p condition number = %8.3e', hcond)); + end + +end diff --git a/Optimization/l1qc_logbarrier.m b/Optimization/l1qc_logbarrier.m new file mode 100644 index 0000000..388529e --- /dev/null +++ b/Optimization/l1qc_logbarrier.m @@ -0,0 +1,116 @@ +% 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, cgres] = 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; + +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 + diff --git a/Optimization/l1qc_newton.m b/Optimization/l1qc_newton.m new file mode 100644 index 0000000..8a25cd2 --- /dev/null +++ b/Optimization/l1qc_newton.m @@ -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 + + + + diff --git a/Optimization/tvdantzig_logbarrier.m b/Optimization/tvdantzig_logbarrier.m new file mode 100644 index 0000000..39c9463 --- /dev/null +++ b/Optimization/tvdantzig_logbarrier.m @@ -0,0 +1,120 @@ +% tvdantzig_logbarrier.m +% +% Solve the total variation Dantzig program +% +% min_x TV(x) subject to ||A'(Ax-b)||_\infty <= epsilon +% +% Recast as the SOCP +% min sum(t) s.t. ||D_{ij}x||_2 <= t, i,j=1,...,n +% <a_{ij},Ax - b> <= epsilon i,j=1,...,n +% and use a log barrier algorithm. +% +% Usage: xp = tvdantzig_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 = tvdantzig_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); +n = round(sqrt(N)); + +% create (sparse) differencing matrices for TV +Dv = spdiags([reshape([-ones(n-1,n); zeros(1,n)],N,1) ... + reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N); +Dh = spdiags([reshape([-ones(n,n-1) zeros(n,1)],N,1) ... + reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N); + +if (largescale) + if (norm(A(x0)-b) > epsilon) + disp('Starting point infeasible; using x0 = At*inv(AAt)*y.'); + AAt = @(z) A(At(z)); + [w, cgres] = 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; +Dhx = Dh*x; Dvx = Dv*x; +t = 1.05*sqrt(Dhx.^2 + Dvx.^2) + .01*max(sqrt(Dhx.^2 + Dvx.^2)); + +% choose initial value of tau so that the duality gap after the first +% step will be about the origial TV +tau = 3*N/sum(sqrt(Dhx.^2+Dvx.^2)); + +lbiter = ceil((log(3*N)-log(lbtol)-log(tau))/log(mu)); +disp(sprintf('Number of log barrier iterations = %d\n', lbiter)); +totaliter = 0; +for ii = 1:lbiter + + [xp, tp, ntiter] = tvdantzig_newton(x, t, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter); + totaliter = totaliter + ntiter; + tvxp = sum(sqrt((Dh*xp).^2 + (Dv*xp).^2)); + + disp(sprintf('\nLog barrier iter = %d, TV = %.3f, functional = %8.3f, tau = %8.3e, total newton iter = %d\n', ... + ii, tvxp, sum(tp), tau, totaliter)); + + x = xp; + t = tp; + + tau = mu*tau; + +end +
\ No newline at end of file diff --git a/Optimization/tvdantzig_newton.m b/Optimization/tvdantzig_newton.m new file mode 100644 index 0000000..68d7148 --- /dev/null +++ b/Optimization/tvdantzig_newton.m @@ -0,0 +1,185 @@ +% tvdantzig_newton.m +% +% Newton iterations for TV Dantzig log-barrier subproblem. +% +% Usage : [xp, tp, niter] = tvdantzig_newton(x0, t0, A, At, b, epsilon, tau, +% newtontol, newtonmaxiter, cgtol, cgmaxiter) +% +% x0,t0 - Nx1 vectors, initial 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. +% +% newtonmaxiter - Maximum number of iterations. +% +% cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix. +% +% cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored +% if A is a matrix. +% +% Written by: Justin Romberg, Caltech +% Email: jrom@acm.caltech.edu +% Created: October 2005 +% + + +function [xp, tp, niter] = tvdantzig_newton(x0, t0, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter) + +largescale = isa(A,'function_handle'); + +alpha = 0.01; +beta = 0.5; + +N = length(x0); +n = round(sqrt(N)); + +% create (sparse) differencing matrices for TV +Dv = spdiags([reshape([-ones(n-1,n); zeros(1,n)],N,1) ... + reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N); +Dh = spdiags([reshape([-ones(n,n-1) zeros(n,1)],N,1) ... + reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N); + +% initial point +x = x0; +t = t0; +if (largescale) + r = A(x) - b; + Atr = At(r); +else + AtA = A'*A; + r = A*x - b; + Atr = A'*r; +end +Dhx = Dh*x; Dvx = Dv*x; +ft = 1/2*(Dhx.^2 + Dvx.^2 - t.^2); +fe1 = Atr - epsilon; +fe2 = -Atr - epsilon; +f = sum(t) - (1/tau)*(sum(log(-ft)) + sum(log(-fe1)) + sum(log(-fe2))); + +niter = 0; +done = 0; +while (~done) + + if (largescale) + ntgx = Dh'*((1./ft).*Dhx) + Dv'*((1./ft).*Dvx) + At(A(1./fe1-1./fe2)); + else + ntgx = Dh'*((1./ft).*Dhx) + Dv'*((1./ft).*Dvx) + AtA*(1./fe1-1./fe2); + end + ntgt = -tau - t./ft; + gradf = -(1/tau)*[ntgx; ntgt]; + + sig22 = 1./ft + (t.^2)./(ft.^2); + sig12 = -t./ft.^2; + sigb = 1./ft.^2 - (sig12.^2)./sig22; + siga = 1./fe1.^2 + 1./fe2.^2; + + w11 = ntgx - Dh'*(Dhx.*(sig12./sig22).*ntgt) - Dv'*(Dvx.*(sig12./sig22).*ntgt); + if (largescale) + h11pfun = @(w) H11p(w, A, At, Dh, Dv, Dhx, Dvx, sigb, ft, siga); + [dx, cgres, cgiter] = cgsolve(h11pfun, w11, cgtol, cgmaxiter, 0); + if (cgres > 1/2) + disp('Cannot solve system. Returning previous iterate. (See Section 4 of notes for more information.)'); + xp = x; tp = t; + return + end + Adx = A(dx); + AtAdx = At(Adx); + else + H11p = Dh'*sparse(diag(-1./ft + sigb.*Dhx.^2))*Dh + ... + Dv'*sparse(diag(-1./ft + sigb.*Dvx.^2))*Dv + ... + Dh'*sparse(diag(sigb.*Dhx.*Dvx))*Dv + ... + Dv'*sparse(diag(sigb.*Dhx.*Dvx))*Dh + ... + AtA*sparse(diag(siga))*AtA; + opts.POSDEF = true; opts.SYM = true; + [dx,hcond] = linsolve(H11p, w11, opts); + if (hcond < 1e-14) + disp('Matrix ill-conditioned. Returning previous iterate. (See Section 4 of notes for more information.)'); + xp = x; tp = t; + return + end + Adx = A*dx; + AtAdx = A'*Adx; + end + Dhdx = Dh*dx; Dvdx = Dv*dx; + dt = (1./sig22).*(ntgt - sig12.*(Dhx.*Dhdx + Dvx.*Dvdx)); + + % minimum step size that stays in the interior + ife1 = find(AtAdx > 0); ife2 = find(-AtAdx > 0); + aqt = Dhdx.^2 + Dvdx.^2 - dt.^2; + bqt = 2*(Dhdx.*Dhx + Dvdx.*Dvx - t.*dt); + cqt = Dhx.^2 + Dvx.^2 - t.^2; + tsols = [(-bqt+sqrt(bqt.^2-4*aqt.*cqt))./(2*aqt); ... + (-bqt-sqrt(bqt.^2-4*aqt.*cqt))./(2*aqt) ]; + indt = find([(bqt.^2 > 4*aqt.*cqt); (bqt.^2 > 4*aqt.*cqt)] & (tsols > 0)); + smax = min(1, min([-fe1(ife1)./AtAdx(ife1); -fe2(ife2)./(-AtAdx(ife2)); tsols(indt)])); + s = (0.99)*smax; + + % backtracking line search + suffdec = 0; + backiter = 0; + while (~suffdec) + xp = x + s*dx; tp = t + s*dt; + rp = r + s*Adx; Atrp = Atr + s*AtAdx; + Dhxp = Dhx + s*Dhdx; Dvxp = Dvx + s*Dvdx; + ftp = 1/2*(Dhxp.^2 + Dvxp.^2 - tp.^2); + fe1p = Atrp - epsilon; + fe2p = -Atrp - epsilon; + fp = sum(tp) - (1/tau)*(sum(log(-ftp)) + sum(log(-fe1p)) + sum(log(-fe2p))); + flin = f + alpha*s*(gradf'*[dx; dt]); + suffdec = (fp <= flin); + s = beta*s; + backiter = backiter + 1; + if (backiter > 32) + disp('Stuck backtracking, returning previous iterate. (See Section 4 of notes for more information.)'); + xp = x; tp = t; + return + end + end + + % set up for next iteration + x = xp; t = tp; + r = rp; Atr = Atrp; + Dvx = Dvxp; Dhx = Dhxp; + ft = ftp; fe1 = fe1p; fe2 = fe2p; f = fp; + + lambda2 = -(gradf'*[dx; dt]); + stepsize = s*norm([dx; dt]); + 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 + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% H11p auxiliary function +function y = H11p(v, A, At, Dh, Dv, Dhx, Dvx, sigb, ft, siga) + +Dhv = Dh*v; +Dvv = Dv*v; + +y = Dh'*((-1./ft + sigb.*Dhx.^2).*Dhv + sigb.*Dhx.*Dvx.*Dvv) + ... + Dv'*((-1./ft + sigb.*Dvx.^2).*Dvv + sigb.*Dhx.*Dvx.*Dhv) + ... + At(A(siga.*At(A(v)))); + + diff --git a/Optimization/tveq_logbarrier.m b/Optimization/tveq_logbarrier.m new file mode 100644 index 0000000..617bf2e --- /dev/null +++ b/Optimization/tveq_logbarrier.m @@ -0,0 +1,118 @@ +% tveq_logbarrier.m +% +% Solve equality constrained TV minimization +% min TV(x) s.t. Ax=b. +% +% Recast as the SOCP +% min sum(t) s.t. ||D_{ij}x||_2 <= t, i,j=1,...,n +% Ax=b +% and use a log barrier algorithm. +% +% Usage: xp = tveq_logbarrier(x0, A, At, b, lbtol, mu, slqtol, slqmaxiter) +% +% 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. +% +% 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. +% +% slqtol - Tolerance for SYMMLQ; ignored if A is a matrix. +% Default = 1e-8. +% +% slqmaxiter - Maximum number of iterations for SYMMLQ; ignored +% if A is a matrix. +% Default = 200. +% +% Written by: Justin Romberg, Caltech +% Email: jrom@acm.caltech.edu +% Created: October 2005 +% + +function xp = tveq_logbarrier(x0, A, At, b, lbtol, mu, slqtol, slqmaxiter) + +largescale = isa(A,'function_handle'); + +if (nargin < 5), lbtol = 1e-3; end +if (nargin < 6), mu = 10; end +if (nargin < 7), slqtol = 1e-8; end +if (nargin < 8), slqmaxiter = 200; end + +newtontol = lbtol; +newtonmaxiter = 50; + +N = length(x0); +n = round(sqrt(N)); + +% create (sparse) differencing matrices for TV +Dv = spdiags([reshape([-ones(n-1,n); zeros(1,n)],N,1) ... + reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N); +Dh = spdiags([reshape([-ones(n,n-1) zeros(n,1)],N,1) ... + reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N); + +% starting point --- make sure that it is feasible +if (largescale) + if (norm(A(x0)-b)/norm(b) > slqtol) + disp('Starting point infeasible; using x0 = At*inv(AAt)*y.'); + AAt = @(z) A(At(z)); + [w,cgres] = cgsolve(AAt, b, slqtol, slqmaxiter, 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)/norm(b) > slqtol) + 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; +Dhx = Dh*x; Dvx = Dv*x; +t = (0.95)*sqrt(Dhx.^2 + Dvx.^2) + (0.1)*max(sqrt(Dhx.^2 + Dvx.^2)); + +% choose initial value of tau so that the duality gap after the first +% step will be about the origial TV +tau = N/sum(sqrt(Dhx.^2+Dvx.^2)); + +lbiter = ceil((log(N)-log(lbtol)-log(tau))/log(mu)); +disp(sprintf('Number of log barrier iterations = %d\n', lbiter)); +totaliter = 0; +for ii = 1:lbiter + + [xp, tp, ntiter] = tveq_newton(x, t, A, At, b, tau, newtontol, newtonmaxiter, slqtol, slqmaxiter); + totaliter = totaliter + ntiter; + + tvxp = sum(sqrt((Dh*xp).^2 + (Dv*xp).^2)); + disp(sprintf('\nLog barrier iter = %d, TV = %.3f, functional = %8.3f, tau = %8.3e, total newton iter = %d\n', ... + ii, tvxp, sum(tp), tau, totaliter)); + + x = xp; + t = tp; + + tau = mu*tau; + +end +
\ No newline at end of file diff --git a/Optimization/tveq_newton.m b/Optimization/tveq_newton.m new file mode 100644 index 0000000..9e71b73 --- /dev/null +++ b/Optimization/tveq_newton.m @@ -0,0 +1,180 @@ +% tveq_newton.m +% +% Newton algorithm for log-barrier subproblems for TV minimization +% with equality constraints. +% +% Usage: +% [xp,tp,niter] = tveq_newton(x0, t0, A, At, b, tau, +% newtontol, newtonmaxiter, slqtol, slqmaxiter) +% +% x0,t0 - 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. +% +% tau - Log barrier parameter. +% +% newtontol - Terminate when the Newton decrement is <= newtontol. +% +% newtonmaxiter - Maximum number of iterations. +% +% slqtol - Tolerance for SYMMLQ; ignored if A is a matrix. +% +% slqmaxiter - Maximum number of iterations for SYMMLQ; ignored +% if A is a matrix. +% +% Written by: Justin Romberg, Caltech +% Email: jrom@acm.caltech.edu +% Created: October 2005 +% + +function [xp, tp, niter] = tveq_newton(x0, t0, A, At, b, tau, newtontol, newtonmaxiter, slqtol, slqmaxiter) + +largescale = isa(A,'function_handle'); + +alpha = 0.01; +beta = 0.5; + +N = length(x0); +n = round(sqrt(N)); +K = length(b); + +% create (sparse) differencing matrices for TV +Dv = spdiags([reshape([-ones(n-1,n); zeros(1,n)],N,1) ... + reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N); +Dh = spdiags([reshape([-ones(n,n-1) zeros(n,1)],N,1) ... + reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N); + +% auxillary matrices for preconditioning +Mdv = spdiags([reshape([ones(n-1,n); zeros(1,n)],N,1) ... + reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N); +Mdh = spdiags([reshape([ones(n,n-1) zeros(n,1)],N,1) ... + reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N); +Mmd = reshape([ones(n-1,n-1) zeros(n-1,1); zeros(1,n)],N,1); + + +% initial point +x = x0; +t = t0; +Dhx = Dh*x; Dvx = Dv*x; +ft = 1/2*(Dhx.^2 + Dvx.^2 - t.^2); +f = sum(t) - (1/tau)*(sum(log(-ft))); + +niter = 0; +done = 0; +while (~done) + + ntgx = Dh'*((1./ft).*Dhx) + Dv'*((1./ft).*Dvx); + ntgt = -tau - t./ft; + gradf = -(1/tau)*[ntgx; ntgt]; + + sig22 = 1./ft + (t.^2)./(ft.^2); + sig12 = -t./ft.^2; + sigb = 1./ft.^2 - (sig12.^2)./sig22; + + w1p = ntgx - Dh'*(Dhx.*(sig12./sig22).*ntgt) - Dv'*(Dvx.*(sig12./sig22).*ntgt); + wp = [w1p; zeros(K,1)]; + if (largescale) + % diagonal of H11p + dg11p = Mdh'*(-1./ft + sigb.*Dhx.^2) + Mdv'*(-1./ft + sigb.*Dvx.^2) + 2*Mmd.*sigb.*Dhx.*Dvx; + afac = max(dg11p); + hpfun = @(z) Hpeval(z, A, At, Dh, Dv, Dhx, Dvx, sigb, ft, afac); + [dxv,slqflag,slqres,slqiter] = symmlq(hpfun, wp, slqtol, slqmaxiter); + if (slqres > 1/2) + disp('Cannot solve system. Returning previous iterate. (See Section 4 of notes for more information.)'); + xp = x; + return + end + else + H11p = Dh'*sparse(diag(-1./ft + sigb.*Dhx.^2))*Dh + ... + Dv'*sparse(diag(-1./ft + sigb.*Dvx.^2))*Dv + ... + Dh'*sparse(diag(sigb.*Dhx.*Dvx))*Dv + ... + Dv'*sparse(diag(sigb.*Dhx.*Dvx))*Dh; + afac = max(diag(H11p)); + Hp = full([H11p afac*A'; afac*A zeros(K)]); + %keyboard + opts.SYM = true; + [dxv, hcond] = linsolve(Hp, wp, opts); + if (hcond < 1e-14) + disp('Matrix ill-conditioned. Returning previous iterate. (See Section 4 of notes for more information.)'); + xp = x; tp = t; + return + end + end + dx = dxv(1:N); + Dhdx = Dh*dx; Dvdx = Dv*dx; + dt = (1./sig22).*(ntgt - sig12.*(Dhx.*Dhdx + Dvx.*Dvdx)); + + % minimum step size that stays in the interior + aqt = Dhdx.^2 + Dvdx.^2 - dt.^2; + bqt = 2*(Dhdx.*Dhx + Dvdx.*Dvx - t.*dt); + cqt = Dhx.^2 + Dvx.^2 - t.^2; + tsols = [(-bqt+sqrt(bqt.^2-4*aqt.*cqt))./(2*aqt); ... + (-bqt-sqrt(bqt.^2-4*aqt.*cqt))./(2*aqt) ]; + indt = find([(bqt.^2 > 4*aqt.*cqt); (bqt.^2 > 4*aqt.*cqt)] & (tsols > 0)); + smax = min(1, min(tsols(indt))); + s = (0.99)*smax; + + % line search + suffdec = 0; + backiter = 0; + while (~suffdec) + xp = x + s*dx; tp = t + s*dt; + Dhxp = Dhx + s*Dhdx; Dvxp = Dvx + s*Dvdx; + ftp = 1/2*(Dhxp.^2 + Dvxp.^2 - tp.^2); + fp = sum(tp) - (1/tau)*(sum(log(-ftp))); + flin = f + alpha*s*(gradf'*[dx; dt]); + suffdec = (fp <= flin); + s = beta*s; + backiter = backiter + 1; + if (backiter > 32) + disp('Stuck backtracking, returning last iterate. (See Section 4 of notes for more information.)'); + xp = x; tp = t; + return + end + end + + % set up for next iteration + x = xp; t = tp; + Dvx = Dvxp; Dhx = Dhxp; + ft = ftp; f = fp; + + lambda2 = -(gradf'*[dx; dt]); + stepsize = s*norm([dx; dt]); + 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(' SYMMLQ Res = %8.3e, SYMMLQ Iter = %d', slqres, slqiter)); + else + disp(sprintf(' H11p condition number = %8.3e', hcond)); + end + +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Implicit application of Hessian +function y = Hpeval(z, A, At, Dh, Dv, Dhx, Dvx, sigb, ft, afac) + +N = length(ft); +K = length(z)-N; +w = z(1:N); +v = z(N+1:N+K); + +Dhw = Dh*w; +Dvw = Dv*w; + +y1 = Dh'*((-1./ft + sigb.*Dhx.^2).*Dhw + sigb.*Dhx.*Dvx.*Dvw) + ... + Dv'*((-1./ft + sigb.*Dvx.^2).*Dvw + sigb.*Dhx.*Dvx.*Dhw) + afac*At(v); +y2 = afac*A(w); + +y = [y1; y2]; diff --git a/Optimization/tvqc_logbarrier.m b/Optimization/tvqc_logbarrier.m new file mode 100644 index 0000000..52cc6b1 --- /dev/null +++ b/Optimization/tvqc_logbarrier.m @@ -0,0 +1,121 @@ +% tvqc_logbarrier.m +% +% Solve quadractically constrained TV minimization +% min TV(x) s.t. ||Ax-b||_2 <= epsilon. +% +% Recast as the SOCP +% min sum(t) s.t. ||D_{ij}x||_2 <= t, i,j=1,...,n +% ||Ax - b||_2 <= epsilon +% and use a log barrier algorithm. +% +% Usage: xp = tvqc_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, tp] = tvqc_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); +n = round(sqrt(N)); + +% create (sparse) differencing matrices for TV +Dv = spdiags([reshape([-ones(n-1,n); zeros(1,n)],N,1) ... + reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N); +Dh = spdiags([reshape([-ones(n,n-1) zeros(n,1)],N,1) ... + reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N); + +% 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, cgres] = 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; +Dhx = Dh*x; Dvx = Dv*x; +t = (0.95)*sqrt(Dhx.^2 + Dvx.^2) + (0.1)*max(sqrt(Dhx.^2 + Dvx.^2)); + +% choose initial value of tau so that the duality gap after the first +% step will be about the origial TV +tau = (N+1)/sum(sqrt(Dhx.^2+Dvx.^2)); + +lbiter = ceil((log((N+1))-log(lbtol)-log(tau))/log(mu)); +disp(sprintf('Number of log barrier iterations = %d\n', lbiter)); +totaliter = 0; +for ii = 1:lbiter + + [xp, tp, ntiter] = tvqc_newton(x, t, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter); + totaliter = totaliter + ntiter; + + tvxp = sum(sqrt((Dh*xp).^2 + (Dv*xp).^2)); + disp(sprintf('\nLog barrier iter = %d, TV = %.3f, functional = %8.3f, tau = %8.3e, total newton iter = %d\n', ... + ii, tvxp, sum(tp), tau, totaliter)); + + x = xp; + t = tp; + + tau = mu*tau; + +end + diff --git a/Optimization/tvqc_newton.m b/Optimization/tvqc_newton.m new file mode 100644 index 0000000..febe8ff --- /dev/null +++ b/Optimization/tvqc_newton.m @@ -0,0 +1,176 @@ +% tvqc_newton.m +% +% Newton algorithm for log-barrier subproblems for TV minimization +% with quadratic constraints. +% +% Usage: +% [xp,tp,niter] = tvqc_newton(x0, t0, A, At, b, epsilon, tau, +% newtontol, newtonmaxiter, cgtol, cgmaxiter) +% +% x0,t0 - 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. +% +% newtonmaxiter - Maximum number of iterations. +% +% cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix. +% +% cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored +% if A is a matrix. +% +% Written by: Justin Romberg, Caltech +% Email: jrom@acm.caltech.edu +% Created: October 2005 +% + +function [xp, tp, niter] = tvqc_newton(x0, t0, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter) + +largescale = isa(A,'function_handle'); + +alpha = 0.01; +beta = 0.5; + +N = length(x0); +n = round(sqrt(N)); + +% create (sparse) differencing matrices for TV +Dv = spdiags([reshape([-ones(n-1,n); zeros(1,n)],N,1) ... + reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N); +Dh = spdiags([reshape([-ones(n,n-1) zeros(n,1)],N,1) ... + reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N); + +if (~largescale), AtA = A'*A; end; + +% initial point +x = x0; +t = t0; +if (largescale), r = A(x) - b; else r = A*x - b; end +Dhx = Dh*x; Dvx = Dv*x; +ft = 1/2*(Dhx.^2 + Dvx.^2 - t.^2); +fe = 1/2*(r'*r - epsilon^2); +f = sum(t) - (1/tau)*(sum(log(-ft)) + log(-fe)); + +niter = 0; +done = 0; +while (~done) + + if (largescale), Atr = At(r); else Atr = A'*r; end + ntgx = Dh'*((1./ft).*Dhx) + Dv'*((1./ft).*Dvx) + 1/fe*Atr; + ntgt = -tau - t./ft; + gradf = -(1/tau)*[ntgx; ntgt]; + + sig22 = 1./ft + (t.^2)./(ft.^2); + sig12 = -t./ft.^2; + sigb = 1./ft.^2 - (sig12.^2)./sig22; + + w1p = ntgx - Dh'*(Dhx.*(sig12./sig22).*ntgt) - Dv'*(Dvx.*(sig12./sig22).*ntgt); + if (largescale) + h11pfun = @(z) H11p(z, A, At, Dh, Dv, Dhx, Dvx, sigb, ft, fe, 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; tp = t; + return + end + Adx = A(dx); + else + H11p = Dh'*sparse(diag(-1./ft + sigb.*Dhx.^2))*Dh + ... + Dv'*sparse(diag(-1./ft + sigb.*Dvx.^2))*Dv + ... + Dh'*sparse(diag(sigb.*Dhx.*Dvx))*Dv + ... + Dv'*sparse(diag(sigb.*Dhx.*Dvx))*Dh - ... + (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; tp = t; + return + end + Adx = A*dx; + end + Dhdx = Dh*dx; Dvdx = Dv*dx; + dt = (1./sig22).*(ntgt - sig12.*(Dhx.*Dhdx + Dvx.*Dvdx)); + + % minimum step size that stays in the interior + aqt = Dhdx.^2 + Dvdx.^2 - dt.^2; + bqt = 2*(Dhdx.*Dhx + Dvdx.*Dvx - t.*dt); + cqt = Dhx.^2 + Dvx.^2 - t.^2; + tsols = [(-bqt+sqrt(bqt.^2-4*aqt.*cqt))./(2*aqt); ... + (-bqt-sqrt(bqt.^2-4*aqt.*cqt))./(2*aqt) ]; + indt = find([(bqt.^2 > 4*aqt.*cqt); (bqt.^2 > 4*aqt.*cqt)] & (tsols > 0)); + aqe = Adx'*Adx; bqe = 2*r'*Adx; cqe = r'*r - epsilon^2; + smax = min(1,min([... + tsols(indt); ... + (-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; tp = t + s*dt; + rp = r + s*Adx; Dhxp = Dhx + s*Dhdx; Dvxp = Dvx + s*Dvdx; + ftp = 1/2*(Dhxp.^2 + Dvxp.^2 - tp.^2); + fep = 1/2*(rp'*rp - epsilon^2); + fp = sum(tp) - (1/tau)*(sum(log(-ftp)) + log(-fep)); + flin = f + alpha*s*(gradf'*[dx; dt]); + 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; tp = t; + return + end + end + + % set up for next iteration + x = xp; t = tp; + r = rp; Dvx = Dvxp; Dhx = Dhxp; + ft = ftp; fe = fep; f = fp; + + lambda2 = -(gradf'*[dx; dt]); + stepsize = s*norm([dx; dt]); + 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 + + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% H11p auxiliary function +function y = H11p(v, A, At, Dh, Dv, Dhx, Dvx, sigb, ft, fe, atr) + +Dhv = Dh*v; +Dvv = Dv*v; + +y = Dh'*((-1./ft + sigb.*Dhx.^2).*Dhv + sigb.*Dhx.*Dvx.*Dvv) + ... + Dv'*((-1./ft + sigb.*Dvx.^2).*Dvv + sigb.*Dhx.*Dvx.*Dhv) - ... + 1/fe*At(A(v)) + 1/fe^2*(atr'*v)*atr; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |