Code covered by the BSD License

# Fractional Matrix Powers with Frechet Derivatives and Condition Number Estimate

### Lijing Lin (view profile)

03 May 2013 (Updated )

Computing matrix power A^p in complex/real arithmetic, with condition number and Frechet derivatives

test_gallery.m
```%TEST_GALLERY
% from GALLERY.
% "The Matrix Function Toolbox" is needed to run these tests;
% it is available from
%  http://www.maths.manchester.ac.uk/~higham/mftoolbox/
clc
clear all
v = ver;
if ~any(strcmp({v(1:end).Name}, 'Matrix Function Toolbox (MFT)'))
['Matrix Function Toolbox (MFT) must be installed.\n' ...
'Obtain it from ' ...
'http://www.maths.manchester.ac.uk/~higham/mftoolbox' ])
end

% p values
pvals = [ 1/52 1/12 1/3 ];
pvals = [pvals -pvals 1-pvals -1+pvals ];
pvals = sort(pvals,'ascend');
np = length(pvals);
j = 0;

for matsize = [10 20]
n = matsize;         % matrix size
A = randn(n);  % test matrix
e = eig(A);
emin = min(e(real(e)==e));
if emin < 0
A = A + eye(n)*(abs(emin)+1);   % avoid negative real e'vals
end
fprintf('\n-------------------- Test on: %g by %g random matrix -------------------\n',n,n);
fprintf('     p       cond     complex algorithm         real algorithm   \n');
fprintf('                          X          F           X         F   \n');

for ii = 1:np  % Outer loop over different p.

p = pvals(ii);
j = j + 1;

E = randn(n);

XX = expm(p*logm(A));  % 'accurate' matrix power.
%         XX = powerm_x(A,p);
XX_big = expm(p*logm([A,E;zeros(n),A]));
%         XX_big = powerm_x([A,E;zeros(n),A],p);
FF = XX_big(1:n,n+1:2*n);  % 'accurate' Frechet derivative.
nrmX = norm(XX,1);
nrmF = norm(FF,1);
normF(j) = nrmF;
normX(j) = nrmX;

eX0(j) = norm(X0-XX,1)/nrmX;
eF0(j) = norm(F0-FF,1)/nrmF;

eX1(j) = norm(X1-XX,1)/nrmX;
eF1(j) = norm(F1-FF,1)/nrmF;

fprintf('%9.1e %9.1e  %9.1e  %9.1e  %9.1e  %9.1e\n' , ...
p,  cond0(j),  eX0(j),  eF0(j), eX1(j), eF1(j))
end  % ii loop over pvals
end

```