Rational Matrix Functions and Rank-One Updates
(Matlab Codes)
Daniel S. Bernstein and Charles F. Van Loan
| Books | Papers | Research | Biographical | Miscellaneous | Home |
| Kry | Computes Krylov Matrix |
| Rational_Update | Computes R(A+uv') from R(A) |
| Pade | Computes specified Pade Approximate of exp(A) |
| TestA | Checks Rational_Update on Pade Examples |
| phi | "Exact" phi vector using expM |
| phi1 | phi vector using powers of a Taylor approximation to exp |
| phi2 | phi-vector using a single Taylor approximatyion to exp |
| TestB | Compares Efficiency and Accuracy of phi1 and phi2 |
| Jac1 | Finite Difference Approximation of phi1's Jacobian |
| Jac2 | Finite Difference Approximation of phi2's Jacobian |
| FastJac2 | Fast Computation of phi2's Jacobian |
| TestC | Compares Jac2 and FastJac2 |
function K = Kry(A,X,j)
% A is an n-by-n matrix, X is an n-by-p matrix, and j is a positive integer.
% K = [ X , A*X , A^2*X , ... , A^j-1*X ]
[n,p] = size(X);
K = X;
for i=2:j
K = [ K A*K(:,(i-2)*p+1:(i-1)*p)];
end
|top|
Update of a Rational Matrix Function
function [XN,YN,XD,YD,X,Y] = Rational_Update(A,u,v,alfa,N,beta,Q,R)
%
% Assume that A is an n-by-n matrix and that u and v are column n-vectors.
% The n-by-n matrices Q and R define a QR factorization of a nonsingular
% matrix D. This matrix and the n-by-n matrix N are polynomials in A:
%
% N = alfa(0)I + alfa(1)A + ... + alfa(r)A^r
% D = beta(0)I + beta(1)A + ... + beta(q)A^q
%
% The column vectors alfa(1:r) and beta(1:q) contain the coefficients,
% with the exception of alfa(0) and beta(0), which are not required. Note
% that alfa and beta can be empty vectors when the underlying degree is 0.
%
% The output matrices XN,YN,XD,YD,X, and Y relate to
%
% N_ = alfa(0)I + alfa(1)(A+uv') + ... + alfa(p)(A+uv')^r
% D_ = beta(0)I + beta(1)(A+uv') + ... + beta(q)(A+uv')^q
%
% in the following way:
%
% N_ = N + XN*YN' XN and YN are n-by-r
% D_ = N + XD*YD' XD and YD are n-by-q
% (D_\N_) = (D\N) + X*Y' X and Y are n-by-d, d = max([r q])
%
% Exceptions occur if r or q is zero:
%
% If r=0, then XN=YN=zeros(n,1).
% If q=0, then XD=YD=zeros(n,1).
% If r=0 and q=0, then X=Y=zeros(n,1).
%
% Abbreviated calls are allowed for the polynomial case (alfa ~= [], beta = []):
%
% [XN,YN] = Rational_Update(A,u,v,alfa)
%
% 4-argument call:
if nargin==4
d = length(alfa);
XN = Kry(A,u,d);
YN = Kry(A'+v*u',v,d)*hankel(alfa);
else
% General call:
% The special case of a constant numerator polynomial:
if alfa==[]
alfa=0;
end
% The special case of a constant denominator polynomial:
if beta==[],
beta=0;
end
[n,n]=size(A);
r=length(alfa);
q=length(beta);
d=max([r q]);
% Compute the Krylov and Hankel matrices.
Kd = Kry(A,u,d);
Ld = Kry(A'+v*u',v,d);
H_alfa = hankel(alfa);
H_beta = hankel(beta);
% Compute the numerator and denominator perturbations:
XN = Kd(:,1:r); YN = Ld(:,1:r)*H_alfa;
XD = Kd(:,1:q); YD = Ld(:,1:q)*H_beta;
% Now get the rational's perturbations:
YN_tilde = [YN zeros(n,d-r)];
YD_tilde = [YD zeros(n,d-q)];
X = R\(Q'*Kd);
M = eye(d,d) + YD_tilde'*X;
Y = (YN_tilde' - (M\YD_tilde') * (R\(Q'*(N + XN*YN'))))';
end
|top|
function [alfa,beta,N,D] = Pade(p,q,A)
%
% Computes the Pade approximation D^{-1}N where
%
% N = N(p,q) = I + alfa(1)A + ... + alfa(p)A^p
%
% D = D(p,q) = I + beta(1)A + ... + beta(q)A^q
%
% A call of the form Pade(p,q) returns only alfa(1:p) and beta(1:q).
%
%
alfa = zeros(p,1);
a = 1;
for k=1:p
a = a*((p-k+1)/k)/(p+q-k+1);
alfa(k) = a;
end
beta = zeros(q,1);
b = 1;
for k=1:q
b = -b*((q-k+1)/k)/(p+q-k+1);
beta(k) = b;
end
if nargin==3
[n,n] = size(A);
N = eye(n,n);
D = eye(n,n);
Apower = A;
for k=1:max([p q])
if k<=p
N = N + alfa(k)*Apower;
end
if k<=q
D = D + beta(k)*Apower;
end
if k<max([p q])
Apower = A*Apower;
end
end
end
|top|
Test A
% TestA.m
% Confirms the correctness of Rational_Update.
% Uses (p,q) Pade approximations to the matrix exponential.
% The A-matrix is randomly chosen and of modest norm so as to
% avoid conditioning.
% Set the size of the problem.
n = 10;
% Set up A and A1 = A + uv':
A = randn(n,n); A = A/(n*norm(A));
u = randn(n,1); u = u/(n*norm(u));
v = randn(n,1); v = v/(n*norm(v));
A1 = A+u*v';
clc
home
disp('Let F = D\N be the (p,q) Pade approximation of exp(A).')
disp('Let F1 = D1\N1 be the (p,q) Pade approximation of
exp(A+uv'').')
disp('N-Error = norm(N1 - (N+XN*YN''))')
disp('D-Error = norm(D1 - (D+XD*YD''))')
disp('F-Error = norm(F1 - (F+X*Y''))')
disp(' ')
disp(' p q N1-Error D1-Error F1-error')
disp('---------------------------------------------')
for p=0:4
for q=0:4
% Generate the (p,q) Pade approximate of A and A1.
[alfa,beta,N,D] = Pade(p,q,A);
[alfa,beta,N1,D1] = Pade(p,q,A1);
[Q,R] = qr(D);
[XN,YN,XD,YD,X,Y] = Rational_Update(A,u,v,alfa,N,beta,Q,R);
% Compute the error in the updated numerator, denominator, and F
eN = norm(N1-(N+XN*YN'));
eD = norm(D1-(D+XD*YD'));
eF = norm((D1\N1)-((D\N)+X*Y'));
disp(sprintf('%2d %2d %10.3e %10.3e %10.3e',p,q,eN,eD,eF))
end
end
disp(' ')
disp(' (Now test 4-argument calls.)')
disp(' ')
q=0;
for p=1:8
% Generate the (p,0) Pade approximate of A and A1.
[alfa,beta,N,D] = Pade(p,q,A);
[alfa,beta,N1,D1] = Pade(p,q,A1);
[XN,YN] = Rational_Update(A,u,v,alfa);
% Compute the error in the updated numerator.
eN = norm(N1-(N+XN*YN'));
eD = 0;
eF = eN;
disp(sprintf('%2d %2d %10.3e %10.3e %10.3e',p,q,eN,eD,eF))
end
Test A Output:
Let F = D\N be the (p,q) Pade approximation of exp(A).
Let F1 = D1\N1 be the (p,q) Pade approximation of exp(A+uv').
N-Error = norm(N1 - (N+XN*YN'))
D-Error = norm(D1 - (D+XD*YD'))
F-Error = norm(F1 - (F+X*Y'))
p q N1-Error D1-Error F1-error
---------------------------------------------
0 0 -0.000e+00 -0.000e+00 -0.000e+00
0 1 -0.000e+00 2.220e-16 6.706e-16
0 2 -0.000e+00 2.231e-16 4.447e-16
0 3 -0.000e+00 2.244e-16 4.465e-16
0 4 -0.000e+00 2.237e-16 5.571e-16
1 0 2.220e-16 -0.000e+00 2.220e-16
1 1 2.220e-16 2.220e-16 6.668e-16
1 2 2.220e-16 2.221e-16 6.667e-16
1 3 2.220e-16 2.231e-16 4.479e-16
1 4 2.220e-16 4.443e-16 5.573e-16
2 0 2.240e-16 -0.000e+00 2.240e-16
2 1 1.120e-16 2.224e-16 6.668e-16
2 2 2.242e-16 2.221e-16 4.461e-16
2 3 1.111e-16 2.221e-16 4.525e-16
2 4 2.223e-16 2.222e-16 6.680e-16
3 0 4.443e-16 -0.000e+00 4.443e-16
3 1 2.222e-16 1.110e-16 4.478e-16
3 2 2.230e-16 2.225e-16 4.540e-16
3 3 2.222e-16 2.232e-16 4.584e-16
3 4 2.221e-16 2.223e-16 6.728e-16
4 0 4.444e-16 -0.000e+00 4.444e-16
4 1 1.479e-17 2.221e-16 6.672e-16
4 2 2.223e-16 2.221e-16 4.528e-16
4 3 2.255e-16 4.441e-16 6.665e-16
4 4 4.441e-16 2.244e-16 2.257e-16
(Now test 4-argument calls.)
1 0 2.220e-16 0.000e+00 2.220e-16
2 0 2.240e-16 0.000e+00 2.240e-16
3 0 4.443e-16 0.000e+00 4.443e-16
4 0 4.444e-16 0.000e+00 4.444e-16
5 0 4.444e-16 0.000e+00 4.444e-16
6 0 6.664e-16 0.000e+00 6.664e-16
7 0 6.664e-16 0.000e+00 6.664e-16
8 0 6.666e-16 0.000e+00 6.666e-16
|top|
function z = phi(A,t,c,b)
%
% Assume that A is n-by-n and that c and b are column n-vectors.
% Assume that t is a length-m vector with 0 < t(1) < ... < t(m).
% Define X(:,0) = b and delk = t(k)-t(k-1) with t(0)=0;
% Computes a column m-vector z
%
% c'*X(:,1)
% z = :
% c'*X(:,m)
%
% where X(:,k) = exp(delk*A)*X(:,k-1).
%
n = length(b); m = length(t);
%
% Generate X.
%
X = zeros(n,m);
for k=1:m
if k==1;
X(:,k) = expm(A*t(1))*b;
else
X(:,k) = expm(A*(t(k)-t(k-1)))*X(:,k-1);
end
end
%
% Form the z-vector from X.
%
Z = c'*X;
z = zeros(m,1);
z = Z(:);
|top|
A Finite Difference Approximation to phi
function z = phi1(A,t,c,b,r)
%
% Assume that A is n-by-n and that c and b are column n-vectors.
% Assume that t is a length-m vector with 0 < t(1) < ... < t(m).
% Define X(:,0) = b and delk = t(k)-t(k-1) with t(0)=0;
% Computes a column m-vector z
%
% c'*X(:,1)
% z = :
% c'*X(:,m)
%
% where X(:,k) = (I + delk*A + ... + (delk*A)^r/r!)*X(:,k-1).
%
n = length(b); m = length(t);
%
% Generate X.
%
X = zeros(n,m); xprev = b; tprev = 0;
for k=1:m
delk = t(k) - tprev;
X(:,k) = xprev;
term = xprev;
for j=1:r
term = (delk/j)*(A*term);
X(:,k) = X(:,k) + term;
end
tprev = t(k);
xprev = X(:,k);
end
%
% Form the z-vector from X.
%
Z = c'*X;
z = zeros(m,1);
z = Z(:);
|top|
Another Finite Difference Approximation to phi
function z = phi2(A,t,c,b,r)
%
% Assume that A is n-by-n and that c and b are column n-vectors.
% Assume that t is a length-m vector with 0 < t(1) < ... < t(m).
% Computes a column m-vector z
%
% c'*X(:,1)
% z = :
% c'*X(:,m)
%
% where X(:,k) = (I + t(k)*A + ... + (t(k)*A)^r/r!)*b.
%
n = length(b);
m = length(t);
%
% Generate X.
%
X = zeros(n,m);
for k=1:m
X(:,k) = b;
term = b;
for j=1:r
term = (t(k)/j)*(A*term);
X(:,k) = X(:,k) + term;
end
end
%
% Form the z-vector from X.
%
Z = c'*X;
z = zeros(m,1);
z = Z(:);
|top|
Test B
% TestB.m
% Compares the Taylor approximations phi1 and phi2 to
% the "true" phi (obtained with "exact" exponentials)
% Set the problem size. (m >= n^2)
m = 10;
n = 3;
A = randn(n,n);
t = sort(rand(m,1));
b = randn(n,1);
c = randn(n,1);
clc
home
disp(sprintf('n = %2d',n))
disp(sprintf('m = %2d',m))
disp(sprintf('norm(A*t(m)) = %6.1f\n\n',norm(A*t(m))));
disp(' r phi1 phi2 phi1 phi2')
disp(' flops flops error error')
disp('---------------------------------------------------')
flops norm(phi-phi2)
for r = 1:10
z0 = phi(A,t,c,b);
flops(0);
z1 = phi1(A,t,c,b,r);
f1 = flops;
flops(0);
z2 = phi2(A,t,c,b,r);
f2 = flops;
e1 = norm(z1-z0);
e2 = norm(z2-z0);
disp(sprintf('%3d %8d %8d %10.3e %10.3e',r,f1,f2,e1,e2))
end
Test B Output:
n = 3
m = 10
norm(A*t(m)) = 1.8
r phi1 phi2 phi1 phi2
flops flops error error
---------------------------------------------------
1 320 310 5.029e-02 1.168e-01
2 570 560 6.387e-03 5.799e-02
3 820 810 5.351e-04 2.087e-02
4 1070 1060 2.839e-05 4.368e-03
5 1320 1310 1.048e-06 6.069e-04
6 1570 1560 2.960e-08 6.231e-05
7 1820 1810 7.210e-10 5.322e-06
8 2070 2060 1.722e-11 4.369e-07
9 2320 2310 4.106e-13 3.744e-08
10 2570 2560 4.480e-15 3.177e-09
|top|
Finite Difference Jacobian of phi1
function J = Jac1(Ac,t,c,b,r,zc)
%
% Finite difference Jacobian of z(A) = phi1(A,t,c,b,r) at A = Ac.
% Ac is n-by-n, c and b are column n-vectors, t is an m-vector,
% r is a positive integer, and zc = z(Ac).
%
% J is an m-by-n^2 matrix. If k = (j-1)n+i with 1<=i<=n and 1<=j<=n,
% then J(:,k) is an approximation of the partial of z(A) with
% respect to A(i,j) at A = Ac.
m = length(t);
n = length(c);
% There are more sophisticated ways to determine the dividied
% difference parameter, but the following will do for now:
delta = sqrt(eps);
J = zeros(m,n^2);
col = 0;
for j=1:n
for i=1:n
col = col+1;
Aij = Ac;
Aij(i,j) = Ac(i,j)+delta;
J(:,col) = (phi1(Aij,t,c,b,r) - zc)/delta;
end
end
|top|
Finite Difference Jacobian of phi2
function J = Jac2(Ac,t,c,b,r,zc)
%
% Finite difference Jacobian of z(A) = phi2(A,t,c,b,r) at A = Ac.
% Ac is n-by-n, c and b are column n-vectors, t is an m-vector,
% r is a positive integer, and zc = z(Ac).
%
% J is an m-by-n^2 matrix. If k = (j-1)n+i with 1<=i<=n and 1<=j<=n,
% then J(:,k) is an approximation of the partial of z(A) with
% respect to A(i,j) at A = Ac.
m = length(t);
n = length(c);
% There are more sophisticated ways to determine the divided
% difference parameter, but the following will do for now:
delta = sqrt(eps);
J = zeros(m,n^2);
col = 0;
for j=1:n
for i=1:n
col = col+1;
Aij = Ac;
Aij(i,j) = Ac(i,j)+delta;
J(:,col) = (phi2(Aij,t,c,b,r) - zc)/delta;
end
end
|top|
function J = FastJac2(A,t,c,b,r)
%
% The Jacobian of z(A) = phi2(A,t,c,b,r) at A.
% A is n-by-n, c and b are column n-vectors, t is an m-vector,
% r is a positive integer.
%
% J is an m-by-n^2 matrix. If k = (j-1)n+i with 1<=i<=n and 1<=j<=n,
% then J(:,k) is the partial of z(A) with
% respect to A(i,j).
m = length(t);
n = length(c);
K = [eye(n,n) A];
L = [eye(n,n) A'];
for j=2:r-1
P = A*K(:,(j-1)*n+1:j*n);
K = [K P];
L = [L P'];
end
c_vec = K'*c;
b_vec = L'*b;
J = zeros(m,n^2);
[alfa,beta] = Pade(r,0);
for k=1:m
% Set up row k
H = hankel(alfa.*(t(k).^((1:r)')));
col = 0;
for j=1:n
z = H*b_vec(j:n:n*r);
for i=1:n
col = col+1;
J(k,col) = (c_vec(i:n:n*r)'*z);
end
end
end
|top|
Test C
% TestC.m
% Jac2 and Jac2Fast
% Set the problem size. (m >= n^2)
n = 6;
m = n^2;
A = randn(n,n);
t = sort(rand(m,1));
b = randn(n,1);
c = randn(n,1);
clc
home
disp(sprintf('n = %2d',n))
disp(sprintf('m = %2d',m))
disp(sprintf('norm(A*t(m)) = %6.1f\n',norm(A*t(m))));
disp('J is computed via Jac2, Jfast is computed vis FastJac2, and')
disp('rho = max(max(abs(J-Jfast)))')
disp(' ')
disp(' r Jac2 FastJac2 rho')
disp(' flops flops ')
disp('---------------------------------------------------')
flops norm(phi-phi2)
for r = 1:n
flops(0)
zc = phi2(A,t,c,b,r);
J = Jac2(A,t,c,b,r,zc);
f = flops;
flops(0);
Jfast = FastJac2(A,t,c,b,r);
f1 = flops;
err = max(max(abs(J-Jfast)));
disp(sprintf('%3d %8d %8d %10.3e %11.3e',r,f,f1,err))
end
Test C Output
n = 6
m = 36
norm(A*t(m)) = 5.7
J is computed via Jac2, Jfast is computed vis FastJac2, and
rho = max(max(abs(J-Jfast)))
r Jac2 FastJac2 rho
flops flops
---------------------------------------------------
1 131870 4762 6.182e-08
2 245090 10314 1.851e-07
3 358310 15798 2.518e-07
4 471530 22146 3.225e-07
5 584750 29358 3.666e-07
6 697970 37434 3.511e-07
|top|