inverse of matrix times vector

Dear all:

The following code computes the inverse of a matrix times a column vector using the LU decomposition. The result is supposed to be a a covariance matrix, i.e., symmetric, positive semi-definite.

% create arrays
D = [-0.8710, 1.1870, -0.4511;
   1, 0, 0
   0, 1, 0];
R = [1; zeros(8,1)];
A = eye(9) - kron(D,D);
% compute inverse of A times vectorized R using LU decomposition
[L, U] = lu(A);
z = L\R;
x = U\z;
% random draw
gamma = reshape(x,3,3);
mvnrnd(zeros(1,3),gamma)

However, MATLAB returns an error that gamma is not symmetric, positive semi-definite. My hunch is that there is a problem with floating point arithmetic when computing the inverse. Other values for D work fine, but this particular example does not.

Any hints will be greatly appreciated.

Thanks, Peter

Answer 1

[Christine Tobler]

The problem is probably that the vector x is not exactly symmetric (because of floating-point errors). You can try just symmetrizing gamma:

 gamma = (gamma + gamma')/2;

Based on what right-hand size is used, I think it's still possible for gamma to be singular.

If you have the Control Toolbox, you could use dlyap instead, which will return a symmetric result, or even dlyapchol to guarantee a positive definite result.

[Peter]

Thanks, Christine. I was hoping for an elegant (i.e., numerically robust) way to solve the problem (maybe the other two commands you suggest are, but I don't know anything about Lyapunov equations). Of course, your solution works like a charm.

[Christine Tobler]

X = dlyap(A,Q) solves the discrete-time Lyapunov equation

 A*X*A' − X + Q = 0

so I think passing in D and reshape(R, [3 3]) should solve your problem.