2 views (last 30 days)

Show older comments

John D'Errico
on 1 Jul 2021

Edited: John D'Errico
on 1 Jul 2021

Um, the solution is not to find Q. Because it is difficult to constrain a matrix to be positive definite. You don't want to try to do that using lsqnonlin (regardless of what they say.)

Instead, you want to find a triangular matrix, I'll call it R, such that R is UPPER triangular. Then compute Q from R.

Q = R' * R

That is, irregardless what the matrix R is, The product, R' * R will ALWAYS be at least positive semi-definite. The semi-definite case comes from when R is itself singular, and you can even prevent that from happening, if you are careful.

So, for example, pick ANY matrix R.

format long g

R = triu(randn(3))

Is Q positive semi-definite?

Q = R'*R

It is by construction symmetric. The eigenvalues are...

eig(Q)

SURPRISE! It works! Actually, the classic test in MATLAB to see if a matrix is SPD, is to apply chol to it. In fact, what we did was we created Q from the cholesky factorization. So chol should simply recover R. Did it?

chol(Q)

Well, it did, but for the signs on some of those elements. That is arbitrary. This is like a sqrt, the sign is indeterminate.

The point is, you need to search for only a upper triangular matrix R that does what you need. You could also solve this for a LOWER triangular matrix. But then, I suppose you need to call it L. :) Inside the objective function, form Q from R. This will be sufficient to span the set of positive definite matrices that you will need to search over. And best of all, this will be eminently solvable using lsqnonlin. You will search over UPPER triangular matrices R. So instead of 9 elements, you need only search over 6 unknowns.

When will there be a problem, and Q is only semi-definite? Again, that happens only when one of the diagonal elements of R is sufficiently small in magnitude compared to the other diagonal elements. You will need to watch for that, if you care that Q must always be strictly positive definite. The case where Q becomes singular from this is equivalent to an ellipsoid that contains zero volume. So effectively, you have an ellipsoid that has compressed into a plane, or worse, into a stright line. Those are the degenerate cases. But your data should make it painfully obvious when that would happen.

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!