Mercurial > matrix-functions

changeset 7:e0817ccb2834 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave. added funm_files/fun_atom.m added sqrtm2.m added toolbox/gecp.m toolbox/see.m: comment wrong call to subplot toolbox/tmtdemo.m: add a cast to double, as eig does not admit bool matrix input removed funm_files/funm_atom.m
author Antonio Pino Robles <data.script93@gmail.com>
date Tue, 26 May 2015 18:14:54 +0200
parents adcd5fc9b9a2
children a587712dcf5f
files funm_files/fun_atom.m funm_files/funm_atom.m sqrtm2.m toolbox/gecp.m toolbox/see.m toolbox/tmtdemo.m
diffstat 6 files changed, 258 insertions(+), 81 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/funm_files/fun_atom.m	Tue May 26 18:14:54 2015 +0200
@@ -0,0 +1,78 @@
+function [F,n_terms] = fun_atom(T,fun,tol,prnt)
+%FUNM_ATOM  Function of triangular matrix with nearly constant diagonal.
+%           [F, N_TERMS] = FUNM_ATOM(T, FUN, TOL, PRNT) evaluates function
+%           FUN at the upper triangular matrix T, where T has nearly constant
+%           diagonal.  A Taylor series is used.
+%           FUN(X,K) must return the K'th derivative of
+%           the function represented by FUN evaluated at the vector X.
+%           TOL is a convergence tolerance for the Taylor series,
+%           defaulting to EPS.
+%           If PRNT ~= 0 trace information is printed.
+%           N_TERMS is the number of terms taken in the Taylor series.
+%           N_TERMS  = -1 signals lack of convergence.
+
+if nargin < 3 | isempty(tol), tol = eps; end
+if nargin < 4, prnt = 0; end
+
+if isequal(fun,@fun_log)   % LOG is special case.
+   [F,n_terms]  = logm_isst(T,prnt);
+   return
+end
+
+itmax = 500;
+
+n = length(T);
+if n == 1, F = feval(fun,T,0); n_terms = 1; return, end
+
+lambda = sum(diag(T))/n;
+F = eye(n)*feval(fun,lambda,0);
+f_deriv_max = zeros(itmax+n-1,1);
+N = T - lambda*eye(n);
+mu = norm( (eye(n)-abs(triu(T,1)))\ones(n,1),inf );
+
+P = N;
+max_d = 1;
+
+for k = 1:itmax
+
+    f = feval(fun,lambda,k);
+    F_old = F;
+    F = F + P*f;
+    rel_diff = norm(F - F_old,inf)/(tol+norm(F_old,inf));
+    if prnt
+        fprintf('%3.0f: coef = %5.0e', k, abs(f)/factorial(k));
+        fprintf('  N^k/k! = %7.1e', norm(P,inf));
+        fprintf('  rel_d = %5.0e',rel_diff);
+        fprintf('  abs_d = %5.0e',norm(F - F_old,inf));
+    end
+    P = P*N/(k+1);
+
+    if rel_diff <= tol
+
+      % Approximate the maximum of derivatives in convex set containing
+      % eigenvalues by maximum of derivatives at eigenvalues.
+      for j = max_d:k+n-1
+          f_deriv_max(j) = norm(feval(fun,diag(T),j),inf);
+      end
+      max_d = k+n;
+      omega = 0;
+      for j = 0:n-1
+          omega = max(omega,f_deriv_max(k+j)/factorial(j));
+      end
+
+      trunc = norm(P,inf)*mu*omega;  % norm(F) moved to RHS to avoid / 0.
+      if prnt
+          fprintf('  [trunc,test] = [%5.0e %5.0e]', ...
+                   trunc, tol*norm(F,inf))
+      end
+      if prnt == 5, trunc = 0; end % Force simple stopping test.
+      if trunc <= tol*norm(F,inf)
+         n_terms = k;
+         if prnt, fprintf('\n'), end, return
+      end
+    end
+
+    if prnt, fprintf('\n'), end
+
+end
+n_terms = -1;
--- a/funm_files/funm_atom.m	Sat May 16 11:39:36 2015 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,78 +0,0 @@
-function [F,n_terms] = fun_atom(T,fun,tol,prnt)
-%FUNM_ATOM  Function of triangular matrix with nearly constant diagonal.
-%           [F, N_TERMS] = FUNM_ATOM(T, FUN, TOL, PRNT) evaluates function
-%           FUN at the upper triangular matrix T, where T has nearly constant
-%           diagonal.  A Taylor series is used.
-%           FUN(X,K) must return the K'th derivative of
-%           the function represented by FUN evaluated at the vector X.
-%           TOL is a convergence tolerance for the Taylor series,
-%           defaulting to EPS.
-%           If PRNT ~= 0 trace information is printed.
-%           N_TERMS is the number of terms taken in the Taylor series.
-%           N_TERMS  = -1 signals lack of convergence.
-
-if nargin < 3 | isempty(tol), tol = eps; end
-if nargin < 4, prnt = 0; end
-
-if isequal(fun,@fun_log)   % LOG is special case.
-   [F,n_terms]  = logm_isst(T,prnt);
-   return
-end
-
-itmax = 500;
-
-n = length(T);
-if n == 1, F = feval(fun,T,0); n_terms = 1; return, end
-
-lambda = sum(diag(T))/n;
-F = eye(n)*feval(fun,lambda,0);
-f_deriv_max = zeros(itmax+n-1,1);
-N = T - lambda*eye(n);
-mu = norm( (eye(n)-abs(triu(T,1)))\ones(n,1),inf );
-
-P = N;
-max_d = 1;
-
-for k = 1:itmax
-
-    f = feval(fun,lambda,k);
-    F_old = F;
-    F = F + P*f;
-    rel_diff = norm(F - F_old,inf)/(tol+norm(F_old,inf));
-    if prnt
-        fprintf('%3.0f: coef = %5.0e', k, abs(f)/factorial(k));
-        fprintf('  N^k/k! = %7.1e', norm(P,inf));
-        fprintf('  rel_d = %5.0e',rel_diff);
-        fprintf('  abs_d = %5.0e',norm(F - F_old,inf));
-    end
-    P = P*N/(k+1);
-
-    if rel_diff <= tol
-
-      % Approximate the maximum of derivatives in convex set containing
-      % eigenvalues by maximum of derivatives at eigenvalues.
-      for j = max_d:k+n-1
-          f_deriv_max(j) = norm(feval(fun,diag(T),j),inf);
-      end
-      max_d = k+n;
-      omega = 0;
-      for j = 0:n-1
-          omega = max(omega,f_deriv_max(k+j)/factorial(j));
-      end
-
-      trunc = norm(P,inf)*mu*omega;  % norm(F) moved to RHS to avoid / 0.
-      if prnt
-          fprintf('  [trunc,test] = [%5.0e %5.0e]', ...
-                   trunc, tol*norm(F,inf))
-      end
-      if prnt == 5, trunc = 0; end % Force simple stopping test.
-      if trunc <= tol*norm(F,inf)
-         n_terms = k;
-         if prnt, fprintf('\n'), end, return
-      end
-    end
-
-    if prnt, fprintf('\n'), end
-
-end
-n_terms = -1;
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/sqrtm2.m	Tue May 26 18:14:54 2015 +0200
@@ -0,0 +1,121 @@
+function [X, alpha, condest] = sqrtm2(A)
+%SQRTM2    Matrix square root.
+%          X = SQRTM2(A) is a square root of the matrix A (A = X*X).
+%          X is the unique square root for which every eigenvalue has
+%          nonnegative real part (the principal square root).
+%          If A is real with a negative real eigenvalue then a complex
+%          result will be produced.
+%          If A is singular then A may not have a square root.
+%          [X, ALPHA, CONDEST] = SQRTM2(A) returns a stability factor ALPHA and
+%          an estimate CONDEST for the matrix square root condition number of X.
+%          The residual NORM(A-X^2,'fro')/NORM(A,'fro') is bounded by
+%          (n+1)*ALPHA*EPS and the Frobenius norm relative error in X is
+%          bounded approximately by N*ALPHA*CONDEST*EPS, where N = MAX(SIZE(A)).
+
+%          References:
+%          N. J. Higham, Computing real square roots of a real
+%             matrix, Linear Algebra and Appl., 88/89 (1987), pp. 405-430.
+%          A. Bjorck and S. Hammarling, A Schur method for the square root of a
+%             matrix, Linear Algebra and Appl., 52/53 (1983), pp. 127-140.
+
+n = max(size(A));
+[Q, T] = schur(A);        % T is real/complex according to A.
+[Q, T] = rsf2csf(Q, T);   % T is now complex Schur form.
+
+% Compute upper triangular square root R of T, a column at a time.
+
+nzeig = length(find(diag(T)==0));
+if nzeig, warning('Matrix is singular and may not have a square root.'), end
+
+R = zeros(n);
+for j=1:n
+    R(j,j) = sqrt(T(j,j));
+    for i=j-1:-1:1
+        s = 0;
+        for k=i+1:j-1
+            s = s + R(i,k)*R(k,j);
+        end
+        if R(i,i) + R(j,j) ~= 0
+           R(i,j) = (T(i,j) - s)/(R(i,i) + R(j,j));
+        end
+    end
+end
+
+if nargout > 1, alpha = norm(R,'fro')^2 / norm(T,'fro'); end
+
+X = Q*R*Q';
+
+if nargout > 2
+
+   if nzeig
+      condest = inf;
+   else
+
+      % Power method to get condition number estimate.
+      warns = warning;
+      warning('off');
+
+      tol = 1e-2;
+      x = ones(n^2,1);    % Starting vector.
+      cnt = 1;
+      e = 1;
+      e0 = 0;
+      while abs(e-e0) > tol*e & cnt <= 6
+         x = x/norm(x);
+         x0 = x;
+         e0 = e;
+         Sx = solve(R, x);
+         x = solve(R, Sx, 'T');
+         e = sqrt(real(x0'*x));  % sqrt of Rayleigh quotient.
+         fprintf('cnt = %2.0f, e = %9.4e\n', cnt, e)
+         cnt = cnt+1;
+      end
+
+      condest = e*norm(A,'fro')/norm(X,'fro');
+      warning(warns);
+
+   end
+
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function x = solve(R, b, tran)
+%SOLVE       Solves Kronecker system.
+%            x = SOLVE(R, b, TRAN) solves
+%                  A*x = b  if TRAN = '',
+%                 A'*x = b  if TRAN = 'T',
+%            where A = KRON(EYE,R) + KRON(TRANSPOSE(R),EYE).
+%            Default: TRAN = ''.
+
+if nargin < 3, tran = ''; end
+
+n = max(size(R));
+x = zeros(n^2,1);
+
+I = eye(n);
+
+if isempty(tran)
+
+   % Forward substitution.
+   for i = 1:n
+       temp = b(n*(i-1)+1:n*i);
+       for j = 1:i-1
+           temp = temp - R(j,i)*x(n*(j-1)+1:n*j);
+       end
+       x(n*(i-1)+1:n*i) = (R + R(i,i)*I) \ temp;
+   end
+
+elseif strcmp(tran,'T')
+
+   % Back substitution.
+   for i = n:-1:1
+       temp = b(n*(i-1)+1:n*i);
+       for j = i+1:n
+           temp = temp - conj(R(i,j))*x(n*(j-1)+1:n*j);
+       end
+       x(n*(i-1)+1:n*i) = (R' + conj(R(i,i))*I) \ temp;
+   end
+
+end
+
+return
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/toolbox/gecp.m	Tue May 26 18:14:54 2015 +0200
@@ -0,0 +1,56 @@
+
+function [L, U, P, Q, rho] = gecp(A)
+%GECP   Gaussian elimination with complete pivoting.
+%       [L, U, P, Q, RHO] = GECP(A) computes the factorization P*A*Q = L*U,
+%       where L is unit lower triangular, U is upper triangular,
+%       and P and Q are permutation matrices.  RHO is the growth factor.
+%       By itself, GECP(A) returns the final reduced matrix from the
+%       elimination containing both L and U.
+
+%       from the Test Matrix Toolbox version 2.0
+%       https://github.com/carandraug/testmatrix/blob/master/Test%20Matrix%20Toolbox/gecp.m
+
+[n, n] = size(A);
+pp = 1:n; qq = 1:n;
+if nargout == 5
+   maxA = norm(A(:), inf);
+   rho = maxA;
+end
+
+for k = 1:n-1
+
+    % Find largest element in remaining square submatrix.
+    % Note: when tie for max, no guarantee which element is chosen.
+    [colmaxima, rowindices] = max( abs(A(k:n, k:n)) );
+    [biggest, colindex] = max(colmaxima);
+    row = rowindices(colindex)+k-1; col = colindex+k-1;
+
+    % Permute largest element into pivot position.
+    A( [k, row], : ) = A( [row, k], : );
+    A( :, [k, col] ) = A( :, [col, k] );
+    pp( [k, row] ) = pp( [row, k] ); qq( [k, col] ) = qq( [col, k] );
+
+    if A(k,k) == 0
+      break
+    end
+
+    A(k+1:n,k) = A(k+1:n,k)/A(k,k);          % Multipliers.
+
+    % Elimination
+    i = k+1:n;
+    A(i,i) = A(i,i) - A(i,k) * A(k,i);
+    if nargout == 5, rho = max( rho, max(max(abs(A(i,i)))) ); end
+
+end
+
+if nargout <= 1
+   L = A;
+   return
+end
+
+L = tril(A,-1) + eye(n);
+U = triu(A);
+
+if nargout >= 3, P = eye(n); P = P(pp,:); end
+if nargout >= 4, Q = eye(n); Q = Q(:,qq); end
+if nargout == 5, rho = rho/maxA; end
--- a/toolbox/see.m	Sat May 16 11:39:36 2015 +0200
+++ b/toolbox/see.m	Tue May 26 18:14:54 2015 +0200
@@ -59,6 +59,6 @@
       end
       text(0,0,'Matrix not square.')
    end
-   subplot;
+%   subplot;
 
 end
--- a/toolbox/tmtdemo.m	Sat May 16 11:39:36 2015 +0200
+++ b/toolbox/tmtdemo.m	Tue May 26 18:14:54 2015 +0200
@@ -217,9 +217,9 @@
 c = []; e = []; j = 1;
 for i=1:matrix(0)
     A = full(matrix(i, 10));
-    if norm(skewpart(A),1)  % If not Hermitian...
+    if norm(skewpart(A),1) 
        c1 = cond(A);
-       eg = eig(A);
+       eg = eig(double(A));
        e1 = max(abs(eg)) / min(abs(eg));
        % Filter out extremely ill-conditioned matrices.
        if c1 <= 1e10, c(j) = c1; e(j) = e1; j = j + 1; end