nikcleju@18: function [xhat, arepr, lagmult] = ArgminOperL2Constrained(y, M, MH, Omega, OmegaH, Lambdahat, xinit, ilagmult, params) nikcleju@18: nikcleju@18: % nikcleju@18: % This function aims to compute nikcleju@18: % xhat = argmin || Omega(Lambdahat, :) * x ||_2 subject to || y - M*x ||_2 <= epsilon. nikcleju@18: % arepr is the analysis representation corresponding to Lambdahat, i.e., nikcleju@18: % arepr = Omega(Lambdahat, :) * xhat. nikcleju@18: % The function also returns the lagrange multiplier in the process used to compute xhat. nikcleju@18: % nikcleju@18: % Inputs: nikcleju@18: % y : observation/measurements of an unknown vector x0. It is equal to M*x0 + noise. nikcleju@18: % M : Measurement matrix nikcleju@18: % MH : M', the conjugate transpose of M nikcleju@18: % Omega : analysis operator nikcleju@18: % OmegaH : Omega', the conjugate transpose of Omega. Also, synthesis operator. nikcleju@18: % Lambdahat : an index set indicating some rows of Omega. nikcleju@18: % xinit : initial estimate that will be used for the conjugate gradient algorithm. nikcleju@18: % ilagmult : initial lagrange multiplier to be used in nikcleju@18: % params : parameters nikcleju@18: % params.noise_level : this corresponds to epsilon above. nikcleju@18: % params.max_inner_iteration : `maximum' number of iterations in conjugate gradient method. nikcleju@18: % params.l2_accurary : the l2 accuracy parameter used in conjugate gradient method nikcleju@18: % params.l2solver : if the value is 'pseudoinverse', then direct matrix computation (not conjugate gradient method) is used. Otherwise, conjugate gradient method is used. nikcleju@18: % nikcleju@18: nikcleju@18: d = length(xinit); nikcleju@18: lagmultmax = 1e5; nikcleju@18: lagmultmin = 1e-4; nikcleju@18: lagmultfactor = 2; nikcleju@18: accuracy_adjustment_exponent = 4/5; nikcleju@18: lagmult = max(min(ilagmult, lagmultmax), lagmultmin); nikcleju@18: was_infeasible = 0; nikcleju@18: was_feasible = 0; nikcleju@18: nikcleju@18: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nikcleju@18: %% Computation done using direct matrix computation from matlab. (no conjugate gradient method.) nikcleju@18: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nikcleju@18: if strcmp(params.l2solver, 'pseudoinverse') nikcleju@18: if strcmp(class(M), 'double') && strcmp(class(Omega), 'double') nikcleju@18: while true nikcleju@18: alpha = sqrt(lagmult); nikcleju@18: xhat = [M; alpha*Omega(Lambdahat,:)]\[y; zeros(length(Lambdahat), 1)]; nikcleju@18: temp = norm(y - M*xhat, 2); nikcleju@18: %disp(['fidelity error=', num2str(temp), ' lagmult=', num2str(lagmult)]); nikcleju@18: if temp <= params.noise_level nikcleju@18: was_feasible = 1; nikcleju@18: if was_infeasible == 1 nikcleju@18: break; nikcleju@18: else nikcleju@18: lagmult = lagmult*lagmultfactor; nikcleju@18: end nikcleju@18: elseif temp > params.noise_level nikcleju@18: was_infeasible = 1; nikcleju@18: if was_feasible == 1 nikcleju@18: xhat = xprev; nikcleju@18: break; nikcleju@18: end nikcleju@18: lagmult = lagmult/lagmultfactor; nikcleju@18: end nikcleju@18: if lagmult < lagmultmin || lagmult > lagmultmax nikcleju@18: break; nikcleju@18: end nikcleju@18: xprev = xhat; nikcleju@18: end nikcleju@18: arepr = Omega(Lambdahat, :) * xhat; nikcleju@18: return; nikcleju@18: end nikcleju@18: end nikcleju@18: nikcleju@18: nikcleju@18: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nikcleju@18: %% Computation using conjugate gradient method. nikcleju@18: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nikcleju@18: if strcmp(class(MH),'function_handle') nikcleju@18: b = MH(y); nikcleju@18: else nikcleju@18: b = MH * y; nikcleju@18: end nikcleju@18: norm_b = norm(b, 2); nikcleju@18: xhat = xinit; nikcleju@18: xprev = xinit; nikcleju@18: residual = TheHermitianMatrix(xhat, M, MH, Omega, OmegaH, Lambdahat, lagmult) - b; nikcleju@18: direction = -residual; nikcleju@18: iter = 0; nikcleju@18: nikcleju@18: while iter < params.max_inner_iteration nikcleju@18: iter = iter + 1; nikcleju@18: alpha = norm(residual,2)^2 / (direction' * TheHermitianMatrix(direction, M, MH, Omega, OmegaH, Lambdahat, lagmult)); nikcleju@18: xhat = xhat + alpha*direction; nikcleju@18: prev_residual = residual; nikcleju@18: residual = TheHermitianMatrix(xhat, M, MH, Omega, OmegaH, Lambdahat, lagmult) - b; nikcleju@18: beta = norm(residual,2)^2 / norm(prev_residual,2)^2; nikcleju@18: direction = -residual + beta*direction; nikcleju@18: nikcleju@18: if norm(residual,2)/norm_b < params.l2_accuracy*(lagmult^(accuracy_adjustment_exponent)) || iter == params.max_inner_iteration nikcleju@18: if strcmp(class(M), 'function_handle') nikcleju@18: temp = norm(y-M(xhat), 2); nikcleju@18: else nikcleju@18: temp = norm(y-M*xhat, 2); nikcleju@18: end nikcleju@18: nikcleju@18: if strcmp(class(Omega), 'function_handle') nikcleju@18: u = Omega(xhat); nikcleju@18: u = sqrt(lagmult)*norm(u(Lambdahat), 2); nikcleju@18: else nikcleju@18: u = sqrt(lagmult)*norm(Omega(Lambdahat,:)*xhat, 2); nikcleju@18: end nikcleju@18: nikcleju@18: %disp(['residual=', num2str(norm(residual,2)), ' norm_b=', num2str(norm_b), ' omegapart=', num2str(u), ' fidelity error=', num2str(temp), ' lagmult=', num2str(lagmult), ' iter=', num2str(iter)]); nikcleju@18: nikcleju@18: if temp <= params.noise_level nikcleju@18: was_feasible = 1; nikcleju@18: if was_infeasible == 1 nikcleju@18: break; nikcleju@18: else nikcleju@18: lagmult = lagmultfactor*lagmult; nikcleju@18: residual = TheHermitianMatrix(xhat, M, MH, Omega, OmegaH, Lambdahat, lagmult) - b; nikcleju@18: direction = -residual; nikcleju@18: iter = 0; nikcleju@18: end nikcleju@18: elseif temp > params.noise_level nikcleju@18: lagmult = lagmult/lagmultfactor; nikcleju@18: if was_feasible == 1 nikcleju@18: xhat = xprev; nikcleju@18: break; nikcleju@18: end nikcleju@18: was_infeasible = 1; nikcleju@18: residual = TheHermitianMatrix(xhat, M, MH, Omega, OmegaH, Lambdahat, lagmult) - b; nikcleju@18: direction = -residual; nikcleju@18: iter = 0; nikcleju@18: end nikcleju@18: if lagmult > lagmultmax || lagmult < lagmultmin nikcleju@18: break; nikcleju@18: end nikcleju@18: xprev = xhat; nikcleju@18: %elseif norm(xprev-xhat)/norm(xhat) < 1e-2 nikcleju@18: % disp(['rel_change=', num2str(norm(xprev-xhat)/norm(xhat))]); nikcleju@18: % if strcmp(class(M), 'function_handle') nikcleju@18: % temp = norm(y-M(xhat), 2); nikcleju@18: % else nikcleju@18: % temp = norm(y-M*xhat, 2); nikcleju@18: % end nikcleju@18: % nikcleju@18: % if temp > 1.2*params.noise_level nikcleju@18: % was_infeasible = 1; nikcleju@18: % lagmult = lagmult/lagmultfactor; nikcleju@18: % xprev = xhat; nikcleju@18: % end nikcleju@18: end nikcleju@18: nikcleju@18: end nikcleju@18: disp(['fidelity_error=', num2str(temp)]); nikcleju@18: if iter == params.max_inner_iteration nikcleju@18: %disp('max_inner_iteration reached. l2_accuracy not achieved.'); nikcleju@18: end nikcleju@18: nikcleju@18: %% nikcleju@18: % Compute analysis representation for xhat nikcleju@18: %% nikcleju@18: if strcmp(class(Omega),'function_handle') nikcleju@18: temp = Omega(xhat); nikcleju@18: arepr = temp(Lambdahat); nikcleju@18: else %% here Omega is assumed to be a matrix nikcleju@18: arepr = Omega(Lambdahat, :) * xhat; nikcleju@18: end nikcleju@18: return; nikcleju@18: nikcleju@18: nikcleju@18: %% nikcleju@18: % This function computes (M'*M + lm*Omega(L,:)'*Omega(L,:)) * x. nikcleju@18: %% nikcleju@18: function w = TheHermitianMatrix(x, M, MH, Omega, OmegaH, L, lm) nikcleju@18: if strcmp(class(M), 'function_handle') nikcleju@18: w = MH(M(x)); nikcleju@18: else %% M and MH are matrices nikcleju@18: w = MH * M * x; nikcleju@18: end nikcleju@18: if strcmp(class(Omega),'function_handle') nikcleju@18: v = Omega(x); nikcleju@18: vt = zeros(size(v)); nikcleju@18: vt(L) = v(L); nikcleju@18: w = w + lm*OmegaH(vt); nikcleju@18: else %% Omega is assumed to be a matrix and OmegaH is its conjugate transpose nikcleju@18: w = w + lm*OmegaH(:, L)*Omega(L, :)*x; nikcleju@18: end