[A_s,b_s,L_p,K,M] = std_form(A,L,b,W)

STD_FORM Transform a general-form reg. problem into one in standard form. [A_s,b_s,L_p,K,M] = std_form(A,L,b) (method 1) [A_s,b_s,L_p,x_0] = std_form(A,L,b,W) (method 2) Transforms a regularization problem in general form min { || A x - b ||^2 + lambda^2 || L x ||^2 } into one in standard form min { || A_s x_s - b_s ||^2 + lambda^2 || x_s ||^2 } . Two methods are available. In both methods, the regularized solution to the original problem can be written as x = L_p*x_s + d where L_p and d depend on the method as follows: method = 1: L_p = pseudoinverse of L, d = K*M*(b - A*L_p*x_s) method = 2: L_p = A-weighted pseudoinverse of L, d = x_0. The transformation from x_s back to x can be carried out by means of the subroutine gen_form.

- nu C:\BrainStorm_2001\PublicToolbox\regutools\nu.m
- pinit C:\BrainStorm_2001\PublicToolbox\regutools\pinit.m

function [A_s,b_s,L_p,K,M] = std_form(A,L,b,W) %STD_FORM Transform a general-form reg. problem into one in standard form. % % [A_s,b_s,L_p,K,M] = std_form(A,L,b) (method 1) % [A_s,b_s,L_p,x_0] = std_form(A,L,b,W) (method 2) % % Transforms a regularization problem in general form % min { || A x - b ||^2 + lambda^2 || L x ||^2 } % into one in standard form % min { || A_s x_s - b_s ||^2 + lambda^2 || x_s ||^2 } . % % Two methods are available. In both methods, the regularized % solution to the original problem can be written as % x = L_p*x_s + d % where L_p and d depend on the method as follows: % method = 1: L_p = pseudoinverse of L, d = K*M*(b - A*L_p*x_s) % method = 2: L_p = A-weighted pseudoinverse of L, d = x_0. % % The transformation from x_s back to x can be carried out by means % of the subroutine gen_form. % References: L. Elden, "Algorithms for regularization of ill- % conditioned least-squares problems", BIT 17 (1977), 134-145. % L. Elden, "A weighted pseudoinverse, generalized singular values, % and constrained lest squares problems", BIT 22 (1982), 487-502. % M. Hanke, "Regularization with differential operators. An itera- % tive approach", J. Numer. Funct. Anal. Optim. 13 (1992), 523-540. % Per Christian Hansen, UNI-C, 05/26/93. % Nargin determines which method. if (nargin==3) % Initialization for method 1. [m,n] = size(A); [p,np] = size(L); if (np~=n), error('A and L must have the same number of columns'), end % Special treatment of the case where L is square. if (p==n) L_p = inv(L); K = []; M = []; A_s = A/L; b_s = b; return end % Compute a QR factorization of L'. [K,R] = qr(full(L')); R = R(1:p,:); % Compute a QR factorization of A*K(:,p+1:n)). [H,T] = qr(A*K(:,p+1:n)); T = T(1:n-p,:); % Compute the transformed quantities. L_p = (R\(K(:,1:p)'))'; K = K(:,p+1:n); M = T\(H(:,1:n-p)'); A_s = H(:,n-p+1:m)'*A*L_p; b_s = H(:,n-p+1:m)'*b; else % Initialization for method 2. [m,n] = size(A); [p,nl] = size(L); nu = n-p; if (nl~=n), error('A and L must have the same number of columns'), end % Special treatment of the case where L is square. if (p==n) L_p = inv(L); A_s = A/L; b_s = b; x_0 = zeros(n,1); K = x_0; % Fix output name. return end % Compute NAA and x_0; [NAA,x_0] = pinit(W,A,b); b_s = b - A*x_0; % Compute the transformed quantities. L1 = inv([[eye(nu),zeros(nu,p)];L]); L1 = full(L1(:,nu+1:n)); L_p = L1 - W*(NAA*L1); A_s = A*L_p; % Fix output name. K = x_0; end

Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Oct 12 12:05:14 2004

Cross-Directory links are: ON