Solve matrix equation with some known values, assuming answer is symmetric positive definite

Question

I'd like to solve $\boldsymbol{[K][U]=[F]}$ for $\boldsymbol{[K]}$ assuming that the answer is symmetric positive definite and contains some known values and that $\boldsymbol{[U]}$ and $\boldsymbol{[F]}$ are known. All the matrices are square and invertible, but because the values in $\boldsymbol{[U]}$ are found by a best fit to a cloud of points the result obtained with K=F*inv(U) is not symmetric positive definite. Does anyone have any suggestions of how to do this, preferably in MatLab? Something like this but for matrices perhaps?

Answer 1

A quite sensible solution (and very commonly used) is simply to solve $K=F U^{-1}$ and symmetrize the result: $\tilde{K}=\frac12(K+K^T)$. This will work well if your $K$ is already close to symmetric. Then you can diagonalize the matrix (guaranteed to work because it's symmetric) and set negative eigenvalues to zero and multiply back. In most cases, these approximations don't do much damage to the result, especially if your data is already full of statistical noise.

The mathematically most correct (but wasteful) way is probably to "unroll" the matrix. You can fill K with unknown variables that already take into account the symmetric property. For instance, [a,b;b,c] for 2x2. Then, write an extended system where K is now a vector of [a,b,c,...] that U acts on. You get a huge system that you can solve with pseudo-inverse (least-squares solution that gets you the components of the matrix that solve the system as exactly as possible). This is probably not what you want but I'm including it so you can see what can be done. It wouldn't be so bad because the action of "U" on the "vector K" can be described as a sparse matrix. I'm not sure but this whole process might return the same as the symmetrization.

Again, the positive-definite property is harder to enforce. The best is to just "round" the negative eigenvalues to zero. The lowest eigenvalues always represent the random noise anyway.

Thirdly, I strongly believe that something is done very clumsily, otherwise you wouldn't have this problem. The symmetric property shouldn't break simply because of bad data: the fitting you describe could possibly take into account the implied symmetry. In most cases, the symmetric property reflects some conservation law, nonorientedness of a graph, or something like this. Is it possible that you can ensure the symmetry before solving this equation? Positiveness is more tricky, but for that case, the only think I can say is to just "kill" the negative eigenvalues.

You could also consider SVD decomposition. There, you get $$K=SWV^T$$ where $W$ is diagonal and positive definite (singular values on the diagonal), and $S$ and $V$ are orthogonal. If you matrix is symmetric and positive definite, you then S=V. This is a very good test actually. $SV^T$ should be close to unit matrix. This is your measure how wrong things are. If you get negative diagonal components, or very strong off-diagonal components, this means that part of the matrix violates your conditions. Also, setting small singular values to zero is guaranteed to be the smallest possible reduction of rank (least damage). If the matrix is well behaved, you could also just set S=V (or the other way around) and symmetrize the matrix this way ($K=SWS^T$ with zeroed small components of $W$).

Anyway, it again all depends on the size of your system, expected rank (are there a lot of negligible modes or is it very stronly specified), and how "wrong" the matrix is.

Answer 2

As pointed about above, your problem is overdetermined, because you are requiring $KU=F$, even though $U$ is not certain. If you are willing to relax this requirement to, say, $\|KU-F\|\leq\epsilon$ for some choice of norm and a positive $\epsilon$, then you have a shot at solving your problem.

If you do something like this, then this problem is a semidefinite program, and is related to a common task of positive semidefinite matrix completion. Any SDP solver, including my Matlab toolbox CVX can solve it. Suppose you have stored the row, column, and value of the known elements of $K$ in vectors i, j, and v, respectively. Then here would be the CVX model:

cvx_begin
    variable K(n,n) semidefinite
    minimize(trace(K))
    norm( K * U - F, 'Fro' ) <= epsilon
    for q = 1 : length(v)
        K(i(q),j(q)) == v(q)
    end
cvx_end

Notice that I've chosen to minimize the trace of $K$, which is the sum of the eigenvalues. In theory you could remove the objective altogether and let the solver CVX uses select whatever value of $K$ it happens to converge to. I've also chosen to use the Frobenius norm to measure the closeness of the $KU\approx F$ fit. Another option is to minimize the misfit:

cvx_begin
    variable K(n,n) semidefinite
    minimize( norm( K * U - F, 'Fro' ) )
    for q = 1 : length(v)
        K(i(q),j(q)) == v(q)
    end
cvx_end