I am working on a stochastic processes project. There, I am trying to solve the Wiener-Hopf equation to find the coefficients $c$ of a linear estimator:
$$Rc = r.$$
Here, $R$ is the autocorrelation matrix of a wide sense stationary process, and it is known that this kind of a matrix is positive definite. So I can always take the inverse and have: $$c = R^{-1}r . $$
But, in MATLAB, $R$ is close to singular when size of $R$ is $30\times 30$, I get error:
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.519038e-18.
So when I simulate the mean square error on the growing size of $n$, it needs to always decrease, but after $n\geq 30$ I have increases as well and sometimes the variance is even negative (due to the imprecise inverse). So I want to fix this issue.
I know that this is a natural error to have due to computational precision and it has been asked a lot so far. But, the matrix has a special form, where $R(i,j)$ is equal to $R_X(|i-j|)$, i.e., the autocovariance between $X_i, X_j$. So the matrix has form:
$$\begin{bmatrix} R_X(0) & R_X(1) & R_X(2) \\ R_X(1) & R_X(0) & R_X(1) \\ R_X(2) & R_X(1) & R_X(0) \\ \end{bmatrix}$$
The diagonal is all ones since $R_X(0)=1$ and the matrix is symmetric. I thought maybe there is a linear-algebra trick since this matrix is a symmetric positive definite matrix with diagonal of all ones. It seems to be a very easy matrix and I suspect if there is an analytic way to take the inverse.
Note1: My process is a real valued process, that's why $R$ is a symmetric matrix. Normally the lower half corresponds to the complex conjugate of the upper half.
Note2: An example: when $n=3$, $R$ looks like: $$\begin{bmatrix} 1.0000 & 0.2906 & -0.4020\\ 0.2906 & 1.0000 & 0.2906 \\ -0.4020 & 0.2906 & 1.0000\end{bmatrix}$$
Note3: I tried to compute the Cholesky lower-half to compute inverse with that, of course, I have the same problem. Then, in order to compute the Cholesky decomposition, I followed the eigenvalue approach, but with a similar reason, I have found that some eigenvalues are negative. I think all these issues are similar in some sense... I also tried to solve a linear system via MATLAB's c = linsolve(R,r) property, but of course it uses the function inv(.) as well and leads to the same issue.
Note4: I think this dsp.se answer is very relevant, but not useful.
Final Note: I am using pinv(.), i.e., generalized Moore-Penrose inverse. It still is not accurate (the error is less weird though).