Rational Matrix Functions and Rank-One Updates

(Matlab Codes)

Daniel S. Bernstein  and Charles F. Van Loan

BooksPapers  |  Research   |  Biographical   | MiscellaneousHome |

Contents

 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

Krylov Matrix Computation

```  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.

[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.

[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|

The Exact Vector phi

```  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|

Fast Jacobian of phi2

```  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);

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|

```

```