I have an expression for a covariance matrix $C$ in terms of the indices $i$ and $j$. In this way I can analytically calculate the elements of my covariance matrix, however when I try to invert $C$ matlab gives a warning about the matrix being close to singular. The inversion therefore doesn't work, by which I mean the multiplying $C$ by the inverse given i do not get the identity. I have tried calculating the pseudo inverse but this also does not work. I have also tried adding a small constant along the diagonal but again the results do not work. In general I'm working with matrices of dimension 1200 but I will give low dimensional example matrix that has the same properties, i.e. the matrices are symmetric about the diagonal and the anti-diagonal
19.9939 19.9954 19.9958 19.9951 19.9933 19.9905 19.9865
19.9954 19.9973 19.9981 19.9978 19.9965 19.9940 19.9905
19.9958 19.9981 19.9993 19.9995 19.9985 19.9965 19.9933
19.9951 19.9978 19.9995 20.0000 19.9995 19.9978 19.9951
19.9933 19.9965 19.9985 19.9995 19.9993 19.9981 19.9958
19.9905 19.9940 19.9965 19.9978 19.9981 19.9973 19.9954
19.9865 19.9905 19.9933 19.9951 19.9958 19.9954 19.9939
As mentioned in the title the matrix isn't positive, however the the negative eigenvalues are very small suggesting that the matrix is not positive only due to machine precision. The negative eigenvalues are $-1.4048e-14$ and $-2.4571e-15$. How can I go about inverting these matrices?
Best Answer
Assuming your matrices are like the illustrated one, which I characterize has having all its coefficients close to a constant array (with the constant $20$), you can employ the Sherman-Morrison Formula. This formula relates the inverse of a matrix to the inverse of a perturbation of that matrix. It has been used in statistics for iterative and online formulas for performing least squares regressions, because the inclusion of each new row of data to the design matrix $\mathbb{X}$ has the effect of adding a rank-one matrix to the "sum of squares and products matrix" $\mathbb{X}^\prime \mathbb{X}$.
To illustrate, the matrix $\mathbb{A}$ in the question can be written in the form
$$\mathbb{A} = \mathbb{Y} + 20 (1\cdot 1^\prime)$$
where the entries in $\mathbb{Y}$ appear to lie between $0.00$ and $-0.01$ and $1$ is the $7$-vector (a column) whose components are all ones. The rank-one matrix $20(1\cdot 1^\prime)$ (all of whose components equal $20$) is the perturbation of $\mathbb{Y}$ that yields $\mathbb{A}$. The formula asserts
$$\mathbb{A}^{-1} = \mathbb{Y}^{-1} - \frac{20}{1 + 20 (1^\prime \mathbb{Y}^{-1} 1)} (\mathbb{Y}^{-1} 1)\cdot(1^\prime \mathbb{Y}^{-1}).$$
This amounts to inverting $\mathbb{Y}$ and adjusting its coefficients in terms of the row (or column) sums of that inverse. Provided $\mathbb{Y}$ can be stably inverted, this might be an attractive approach.
As I mentioned in a comment, it is probably unnecessary to invert $\mathbb{A}$ at all. For Fisher scoring, for instance, where the neat mathematical formulas have expressions of the form "$x = \mathbb{A}^{-1}b$" (where typically $b$ is the value of the score function), you should be solving the system
$$\mathbb{A}x = b$$
instead of pre-inverting $\mathbb{A}$ and multiplying its inverse by $b$. This tends to lead to more stable and accurate results.