Adding Elements to Diagonal of Symmetric Matrix to Ensure Positive Definiteness.

Question

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}$.

Answer 1

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

Answer 2

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