| [ KAIST ] in KIDS 글 쓴 이(By): Gatsbi (궁금이) 날 짜 (Date): 2000년 8월 4일 금요일 오후 04시 53분 47초 제 목(Title): Re: [Q] matrix square root 구하는 법 sqrtm은 singular value decomposition (SVD)을 이용한 것입니다. Numerical Recipes in C를 참조하시면 쏘스코드까지 제공되어 있습니다. 수학적인 내용은 선형대수 책에 나와 있습니다. � type sqrtm function [S,esterr] = sqrtm(A) %SQRTM Matrix square root. % S = SQRTM(A) is the matrix square root of A. % Complex results are produced if A has negative eigenvalues. % A warning message is printed if the computed S*S is not close to A. % % [S,esterr] = sqrtm(A) does not print any warning message, but % returns an estimate of the relative residual, norm(S*S-A)/norm(A). % % If A is real, symmetric and positive definite, or complex, Hermitian % and positive definite, then so is the computed matrix square root. % % Some matrices, like A = [0 1; 0 0], do not have any square roots, % real or complex, and SQRTM cannot be expected to produce one. % % See also EXPM, LOGM, FUNM. % CBM, 7-12-92. % Copyright (c) 1984-98 by The MathWorks, Inc. % $Revision: 5.5 $ $Date: 1997/12/05 03:31:41 $ % First try Parlett's method directly. [S,esterr] = funm(A,'sqrt'); tol = 1000*eps; % If funm's error estimate is small, accept the result. if esterr >= tol % Use norm of residual instead of funm's crude error estimate. esterr = norm(S*S-A,1)/norm(A,1); if esterr >= tol | ~isfinite(esterr) [n,n] = size(A); I = eye(n,n); % Try again with a not-quite-random rotation. [J,K] = meshgrid(1:n); R = orth(I + J + K); [S,ignore] = funm(R*A*R','sqrt'); S = R'*S*R; if norm(imag(S),1) <= 1000*tol*norm(S,1), S = real(S); end % One step of improvement. S = (S + A/S)/2; esterr = norm(S*S-A,1)/norm(A,1); end end if (esterr >= tol | ~isfinite(esterr)) & (nargout < 2) warning(['SQRTM appears inaccurate. esterr = ',num2str(esterr)]) end |