5
$\begingroup$

Questions:
1. Is there a simple, interpretable way to determine the distance/closeness of a matrix to being not positive (semi-)definite?

2. Alternatively: how can I systematically create matrices that are just barely positive (semi-)definite?

Background: I've been studying the performance of different estimators used in structural equation modeling (like ML and likes). My focus now is on their convergence behavior when these estimators are based on variance-covariance (VCV) matrices that are close to being invalid VCV's (= not positive (semi-)definite). Up until now, I created them rather crudely by trial and error.

$\endgroup$
2
  • 2
    $\begingroup$ Why not use the inverse of the condition number of the matrix? $\endgroup$ Mar 23, 2017 at 22:01
  • 1
    $\begingroup$ @RodrigodeAzevedo: +1 That would work but it would be slower than a Cholesky and it would be more opaque than using the whole spectrum as both me and terhorst suggest. Also it would potentially be system-dependant as some implementations or samples might break down at different condition numbers; looking at the eigenvalues is more obvious where you should get a break-point (while with c.n. is it $10^{15}$? $10^{17}$? Depends!) You will be just looking to the ratio of $frac{\lambda_{max}}{\lambda_{min}}$ (which is a good idea). $\endgroup$
    – usεr11852
    Mar 23, 2017 at 22:31

2 Answers 2

3
$\begingroup$

As Wikipedia says: "A Hermitian $n \times n$ matrix $A$ is defined as being positive-definite (PD) if the scalar $b^T A b$ is positive for every non-zero column vector $b$ of $n$ real numbers". In addition, $A$ can be equivalently be defined as PD in terms of its eigenvalues $\lambda$ (all being positive) and in terms of its Cholesky decomposition $A = LL^T$ (existing).

  1. The fastest way to check therefore if a given matrix $A$ is PD is to check if $A$ has a Cholesky decomposition. If Cholesky decomposition fails, then $A$ is not a PD matrix. Given that $A$ is PD we expect all the diagonal elements of $L$ to be real and strictly positive. The closer they are to $0$ the closer the matrix $A$ is to not being PD. (This relates closely to a fourth characteristics of PD matrices which is to have positive leading principal minors.) More stringently you can also check the eigenvalues of it and see if they are all positive. The closer they are to $0$ the closer $A$ is to not being PD. This also allows a direct interpretation of "why $A$ is not PD". As an eigenvalue gets closer and closer to $0$ this means that the column space spanned by $A$ is shrinking and at some point one of its dimensions collapses onto another, ie. we have at least one column that is the linear combination of other columns in matrix $A$.

  2. Given a matrix $A$ is PD the fastest way to get a matrix "close to" it that it is not PD will be to set on of the diagonal elements of $L$ to 0 (or something very small). That's one way but it is a bit blunt. You will effectively invalidate one or more of $A$'s leading principal minors. A second way that is more carefully examined would be to use the singular value decomposition $A = USV^T$ and set one of the diagonal elements holding the singular values in $S$ to $0$. This will guarantee that the reconstructed matrix $A^{s_n=0}$ will be the closest matrix to $A$ in term of Frobenius norm (given you choose to "invalidate" the smallest singular value $s$). This method will directly relate to characterisation of $A$ as being PD based on its eigenvalues; in the case of a covariance matrix $A$ the singular values of $A$ equate its eigenvalues so when we set one of the singular values to $0$ we effectively set one of the eigenvalues to $0$. Addendum: We can also use this eigenvalue-based methodology to directly to create near PD matrices without a starting PD matrix $A$ (which is actually the general case). Given any $n \times p$ matrix $B$ we calculate the singular value decomposition $B = USV^T$, set one or more of the diagonal values of $S$ to (near) $0$ values and calculate $B^{s_n = 0} = V S V^T$; $B^{s_n = 0}$ will be PSD.

In R that would be something like:

n = 1000;
p = 10;
set.seed(1);
B = matrix(rnorm(n * p), ncol=p);
svdB = svd(B);
B0 = svdB$v %*% diag(c(svdB$d[1:p-1], .Machine$double.eps)) %*% t(svdB$v) 

Notice that floating-point arithmetic issues will kick-in so even if you set one of the diagonal element of $L$ or $S$ to $0.0000000001$ or something very small, you might still end up with non PD matrix. This will relate to the magnitude of the other elements in $L$ or $S$ (and probably in the numerical linear algebra used (OpenBLAS, MKL, etc.)).

$\endgroup$
3
$\begingroup$

For a real symmetric matrix $M=P\,\text{diag}(\lambda_1,\dots,\lambda_n)\, P^T$, the nearest (in Frobenius norm) PSD matrix is obtained by thresholding the eigenvalues at 0:$$\text{Proj}_{\mathcal{S}_n^+}(M)=P\,\text{diag}\big(\max(\lambda_1,0),\dots,\max(\lambda_n,0)\big)\, P^T.$$ See e.g. here. Hence, one interpretable measure of $M$'s distance from being PSD is simply the distance of its spectrum from the positive orthant.

So, to create almost PSD matrices, take any orthonormal matrix $P$ and then set $M=PDP^T$ for $D=\text{diag}(\lambda_1,\dots,\lambda_n)$ and the $\lambda_i$ are near to the edge of said orthant. For example, in R:

n = 10
P = qr.Q(qr(matrix(rnorm(n * n), nrow=n)))
D = diag(rnorm(n, sd=.1)^2)
M = P %*% D %*% t(P)
$\endgroup$
1
  • 1
    $\begingroup$ +1, this is a good answer. I assumed that the OP already had a covariance matrix and wanted to push it being not PD. Just two comments on your code: 1. Use set.seed to ensure reproducibility. 2. The matrix M is far from "almost PSD" in almost all cases. PSD usually appeaars due to numerical precision issues. You want to "magnify" these issues by having really small eigenvalues as well as large condition numbers. You want something like: D = diag(c(rep(10^8,9),0.1^7)) to ensure you are "almost PSD". The current D is far too lenient. $\endgroup$
    – usεr11852
    Mar 23, 2017 at 21:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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