PINV_SOL - Compute pseudoinverse solution. x = pinv_sol(A,b,method) Returns the pseudoinvers solution (= the minimum norm solution) to the linear equation Ax = b.
0001 function x = pinv_sol(A,b,method) 0002 0003 % PINV_SOL - Compute pseudoinverse solution. 0004 % 0005 % x = pinv_sol(A,b,method) 0006 % 0007 % Returns the pseudoinvers solution (= the minimum norm solution) to 0008 % the linear equation Ax = b. 0009 0010 if nargin < 3 0011 method = 1; 0012 end 0013 0014 switch(method) 0015 case 1 0016 % ----------------------- 0017 % Use QR decomposition. 0018 % ----------------------- 0019 if isempty(A) 0020 x = []; 0021 return 0022 end 0023 0024 [Q1,R1,E] = qr(A',0); 0025 % Account for rank deficiency. The tolerance is 10 x tolerance 0026 % used in \. See qr documentation. 0027 tol = max(size(A))*10*eps*abs(R1(1,1)); 0028 % Compute numerical rank. 0029 r = sum(abs(diag(R1))>tol); 0030 % Compute pseudoinverse solution. 0031 x = Q1(:,1:r)*(R1(1:r,:)'\b(E)); 0032 case 2 0033 % ------------------------------------- 0034 % Use MATLAB's pseudoinverse command. 0035 % ------------------------------------- 0036 x = pinv(A) * b; 0037 case 3 0038 % ------------------------ 0039 % Use SVD decomposition. 0040 % ------------------------ 0041 % Adopted from pinv.m. 0042 [n,m] = size(A); 0043 if n > m 0044 % Ax = b is overdetermined. 0045 [U,S,V] = svd(A,0); 0046 else 0047 % Ax = b is underdetermined. 0048 [V,S,U] = svd(A',0); 0049 end 0050 s = diag(S); 0051 r = sum(s > .0000001); 0052 Sr_pinv = diag(1./s(1:r)); 0053 x = V(:,1:r)*Sr_pinv*U(:,1:r)'*b; 0054 end