inverse of matrix times vector
2 views (last 30 days)
Show older comments
Peter
on 21 Mar 2018
Commented: Christine Tobler
on 22 Mar 2018
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
0 Comments
Accepted Answer
Christine Tobler
on 21 Mar 2018
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.
2 Comments
Christine Tobler
on 22 Mar 2018
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.
More Answers (0)
See Also
Categories
Find more on Matrix Computations in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!