Need matrix inverse for IAA estimation

Hi there,

I know that it is usually not correct to take the inverse of a matrix (\ is much better) but I DO NOT have the usual equation Ax = B. Instead I have a matrix R, and I need

alpha(i,j) = (y(i,j)'*R^-1*d) / (y(i,j)'*R^-1*y(i,j))

y(i,j) is 192*1, R is 192*192, d is 192*1, alpha(i,i) is 1*1

where R = sumi(sumj(|alpha(i,j)|^2*(y(i,j)*Y(i,j)'))) -> this matrix is 192*192

This is iterative for IAA estimation. (find alpha, find R, find alpha, etc...)

Is there a good way to get R^-1 or alpha directly from R? inv just gives me garbage, since it is returning nearly singular matrices, and pinv gives me something that seems close to correct, but the amplitude of alpha is bonkers.

By the way, R is Toeplitz, positive-semidefinite.

Thanks

Answer 1

[Matt J]

Presumably, R is supposed to be initialized at something distinctly non-singular. Or else alpha_ij are supposed to induce non-singularity. Are the y_ij all linearly independent? The Y_ij as well? If so, you should be initializing all the alpha_ij strictly non-zero. If not, I suspect you have bad input data.

[Joseph]

Hi Matt,

Thanks for the reply. I am simulating a Radar in MatLab with several different reconstruction methods (same radar model otherwise, so same data). MMSE and Compressed Sensing recover the image with no problem, so the data is correct. Yes, y(i,j) are all independent since I am transmitting random waveforms.

Thanks, Joseph

[Matt J]

Then I guess I've answered your question? As I've already said, if you start with all alpha_ij>0 and all your Y_ij and y_ij data are independent, the corresponding R should be non-singular, at least at the initial iterations.

If your R converge to something singular at higher iterations, then you need to better understand this algorithm and why it works (we obviously cannot). I can see that, for example, if d=0, you will converge immediately to alpha_ij=0 and R=0 which is singular. Presumably, though, d is supposed to be chosen in a smarter way and the algorithm is supposed to preserve the non-singularity of R at all iterations. If that's the expected property, you can implement the update formula with backslash

alpha_ij = (yij'*(R\d)) / (y_ij'*(R\y_ij))

[Joseph]

The definition of the first alpha is alpha_ij = y_ij'*d/(y_ij'*y_ij) Afterward, R and alpha are calculated using the formula in the first post. The backslash method gives the same solution as using pinv of R, but takes MUCH longer, since pinv (R) only needs to be calculated once for all i,j, but the backslash needs to be calculated every time.

[Matt J]

No, the faster thing (than PINV) would be to concatenate all the y_ij into one matrix Ymat, with each y_ij forming one column. Then you would do

 Q=R\Ymat

Each column of Q would be R\y_ij.

For that matter, you also shouldn't be calculating all the quadratic forms y_ij*R^-1*y_ij individually for each ij. Instead, you should be doing

 quadforms = sum((R\Ymat).*Ymat,1);

This will give you all the quadratic forms at once.

The definition of the first alpha is alpha_ij = y_ij'*d/(y_ij'*y_ij)

Hopefully, you understand now why that may or may not work.