8
$\begingroup$

I have a symmetric matrix $A$, which has zeroes all along the diagonal i.e. $A_{ii}=0$.

I cannot change the off diagonal elements of this matrix, I can only change the diagonal elements. I need this matrix to be positive definite.

One way to ensure this is as follows:

Let $\lambda'$ by the absolute value of the most negative eigenvalue and transform $A\mapsto A + \lambda'I_{na}$. Notice this leaves the off-diagonal elements unchanged, but now it is positive definite:

$(A+\lambda'I_{na})x_{i}=(\lambda_{i}+\lambda')x_{i}=\lambda_{i}^{new}x_{i}$,

where $(x_{i},\lambda_{i})$ are an eigenpair. The eigenvalues of the new matrix formed by adding $\lambda'$ to the diagonal are all positive.

I fear that this solution is sub-optimal - in my application I want to add only as much as I need to and no more. I would like to know if I can ensure positive definiteness by more generally performing

$A+D$

where $D$ is some diagonal matrix.

[Application: Statistics. $A$ is related to the covariance of some augmented data. I want it to be as small as possible so as to reduce the leverage of the missing data.

EDIT: Changed $X_{a}^{T}X_{a}$ to $A$. I was ahead of myself, once I get my desired positive definite matrix I want to set $A_{pd}=A+D=X^{T}_{a}X_{a}$ and take the cholesky decomposition to get $X_{a}$.

$\endgroup$
  • 3
    $\begingroup$ The smallest not necessarily diagonal matrix you can add to make the matrix $X_a^TX_a$ positive semidefinite is $S=\sum -\lambda_ix_ix_i^T$ where $(x_i,\lambda_i)$ range over only the negative eigenpairs of $X_a^TX_a$. What you are looking for is a diagonal matrix $D$ that "dominates" $S$, in the sense that $D-S$ is positive semidefinite. Not sure if there is a closed-form solution, but I believe it can be solved numerically using semidefinite programming. $\endgroup$ – user856 Feb 5 '14 at 19:57
  • 3
    $\begingroup$ Any matrix of the form $X^TX$ is already positive semidefinite. Furthermore, if it has a zero on the diagonal, that the whole row and column of that zero can contain only zeroes. Hence, your matrix $X_a^TX_a$ is a zero matrix. Assuming this was not your intention, please amend the question, so we can provide a proper answer. $\endgroup$ – Vedran Šego Feb 5 '14 at 20:28
  • $\begingroup$ You are right, please ignore $X^{T}X$ and just consider it a matrix $A$. Please see my edit in main post. $\endgroup$ – Lindon Feb 5 '14 at 20:55
  • 1
    $\begingroup$ Please, can you also clarify what do you mean by "I want it to be as small as possible"? Is it in terms of some norm (in $2$-norm, for example, the $\lambda$-shift you have described is optimal; see N. J. Higham. "Computing the polar decomposition - with applications", SIAM J.Sci. Statist. Comput., 7(4):1160–1174, Oct. 1986.) or how else do you define this? $\endgroup$ – Vedran Šego Feb 5 '14 at 23:17
  • 1
    $\begingroup$ @Lindon Your shift will produce a singular semidefinite matrix, so $\log \det(A+D) = \log 0 = -\infty$. Can't get smaller than that. $\endgroup$ – Vedran Šego Feb 13 '14 at 18:33
4
$\begingroup$

Yes, as Rahul stated in the comments, this is a semidefinite program, and a relatively straightforward one at that. In fact, it's very similar to the so-called educational testing problem: $$\begin{array}{ll}\text{maximize} & \textstyle\sum_i D_{ii} \\ \text{subject to} & \Sigma - D \succeq 0 \\ & D \succeq 0 \end{array}$$ In the ETP, $\Sigma$ is already positive semidefinite, and you're subtracting as large of a diagonal matrix as possible without changing that. In contrast, your problem is $$\begin{array}{ll}\text{minimize} & \textstyle\sum_i D_{ii} \\ \text{subject to} & \Sigma + D \succeq 0 \\ & D \succeq 0\end{array}$$ But not surprisingly the methods for solving these problems are extremely similar. Of course, this assumes that you like $\sum_i D_{ii}$ as a measure of how much you are perturbing the matrix. You could also consider $\max_i D_{ii}$, $\sum_i D_{ii}^2$, etc.; as long as the measure remains convex, the problem is readily solved.

Any solver supporting semidefinite programming can handle this. Some software I wrote for MATLAB, CVX, makes this easy; so would similar software YALMIP (also for MATLAB) and CVXOPT for Python. In CVX, your model looks like this:

n = size(A,1);
cvx_begin sdp
    variable d(n)
    minimize(sum(d))
    subject to
        A + diag(d) >= 0
        d >= 0
cvx_end
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Cheers Michael (+1 btw). A first look at SDP is quite overwhelming and I was trying to find an example very close to my own, so introducing me to the educational testing problem was a big help. I'm going to try your code out and will let you know the results :) $\endgroup$ – Lindon Feb 6 '14 at 17:31
  • $\begingroup$ Ah, Michael, an interesting result. For my 8x8 matrix there is one negative eigenvalue and that is -922.8676. I ran your code and the output (rounded) I get for D is d =(0, 919,918, 924, 922, 918, 926, 934), which, except for 0, is pretty much what I was originally doing by adding the identity multiplied by the absolute value of the most negative eigenvalue. I guess I don't know where the 0 is coming from, that requires some thought. It seems Vedran Šego 's response that the lambda-shit is optimal is consistent with the output of your code (except the 0). $\endgroup$ – Lindon Feb 6 '14 at 18:00
  • 1
    $\begingroup$ Any chance you could just offer up your 8x8 matrix values here somewhere? I would be happy to run CVX here on it and see what the deal is. $\endgroup$ – Michael Grant Feb 6 '14 at 20:23
  • $\begingroup$ Hey Michael, cheers, I've uploaded a matlab matrix file which can be found here lindonslog.com/example_code/Amatrix I haven't tried it out yet on a matrix that has more than one negative eigenvalue. That result I still have yet to look at. Also, I'm wondering why you have that d>=0 as a constraint? All I care about is that A+D being positive definite $\endgroup$ – Lindon Feb 6 '14 at 23:32
  • 1
    $\begingroup$ Can't be done with CVX, as that's not a convex model. Maximizations require concave objectives. $\endgroup$ – Michael Grant Oct 23 '14 at 3:19
2
$\begingroup$

This is one minimal adjustment to $A$ to make it positive definite, and you get the $LDL^\top$ decomposition in the process:

In computing $L$ and $D$ at the step where you calculate $D_{jj}$, if that value is too small (smaller than some pre-selected $\epsilon>0$), set it to $\epsilon$ and continue.

You need only add just enough to to ensure $D_{jj}$ is positive at each step.

Now $LDL^\top=A+E$ where $E$ is a diagonal matrix with positive entries. If you pre-pivot your matrix $A$ to put the least diagonally dominant rows at the bottom (and columns to the right) you may get better results than without that preconditioning.

Here is some crude MATLAB code that does the trick:

function [L,D] = modifiedLDLT(A,epsilon)
    % http://mathforcollege.com/nm/mws/gen/04sle/mws_gen_sle_txt_cholesky.pdf

    n = size(A,1);

    % check for valid (square, symmetric) input
    m = size(A,2);
    assert(m==n);
    for i=1:n
        for j=i:n
            assert(A(i,j)==A(j,i), 'not symmetric');
        end
    end
    assert(epsilon>0);

    L = zeros(n,n);
    D = zeros(n,n);

    for j=1:n
        Lsum = 0;
        for k=1:(j-1)
            Lsum = Lsum + L(j,k)*L(j,k)*D(k,k);
        end
        D(j,j) = A(j,j)-Lsum;
        if D(j,j) < epsilon
            D(j,j) = epsilon;
        end

        L(j,j) = 1;
        for i=(j+1):n
            Lsum = 0;
            for k=1:(j-1)
                Lsum = Lsum + L(i,k)*D(k,k)*L(j,k);
            end
            L(i,j) = 1/(D(j,j)) * (A(i,j)-Lsum);
        end
    end
end %function
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ This is great! No sense in solving an SDP (which would do tens of Choleskys) unless you're not OK with this particular measure of closeness. $\endgroup$ – Michael Grant Feb 6 '14 at 4:04
  • $\begingroup$ It's a nice idea Jeff, I coded it up but it doesn't work. At least, not as far as I can tell. For small epsilon it is not stable, and by small I mean epsilon=0.1 or even epsilon=10. Elements get really large and LL* differs on the off-diagonal elements from A, which fails my criteria. For epsilon=10, the off-diagonal elements of LL* and A are the same, but LL* has a huge diagonal, in which case my initial solution of adding the Identity times the absolute value of the most negative eigenvalue is much better. I've provided some code: lindonslog.com/example_code/stackexchange.r $\endgroup$ – Lindon Feb 6 '14 at 17:28
  • $\begingroup$ OK, well, if it doesn't work, then it's not as great as I thought :-) I am familiar with "self-correcting" Cholesky code that provides a sort of "incomplete factorization" for matrices that are just barely indefinite, or appear that way due to roundoff errors in the Cholesky. In fact such code can be very useful for solving semidefinite programs. But perhaps the approach fails if you need to make large modifications. $\endgroup$ – Michael Grant Feb 6 '14 at 20:26
  • 1
    $\begingroup$ I'm unsure why it wouldn't work. Admittedly when I used it I wasn't concerned with keeping the off-diagonals unchanged, but it did give a distinctly diagonal $D$. I will write some code of my own to check that it works, or confirm that it doesn't. $\endgroup$ – Jeff Snider Feb 7 '14 at 2:43
  • $\begingroup$ I have confirmed that while it "works" in the sense of giving a diagonal $D$, the values of $D$ seem to grow exponentially as you proceed down the diagonal, and it ends up being much larger than simply adding $-\lambda_i I$ for the most negative eigenvalue $\lambda_i$. Now I will dig into my old records and see if I can recover precisely what I was doing previously. $\endgroup$ – Jeff Snider Feb 8 '14 at 20:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.