# Iteratively solving a sparse, ill-conditioned system

I have a sparse (density = 0.2%), ill-conditioned system that I am trying to solve, with no luck.

Background

I have a sequence of sampled data, where two of every 8 samples have been zeroed due to a bug. The sequence is two sinusoids in the presence of white noise.

What I'd like to do is determine the Signal to Noise Ratio of the original sequence. I understand that there is information that is forever lost, and so I won't get an exact answer, but that's acceptable, as long as I can get something reasonably close using the existing data, without any filtering (i.e. no interpolation).

Strategy

Having a sequence with missing samples as I do is like multiplying the original sequence with a square wave (75% duty cycle, in my case). In the frequency domain, this is equivalent to a periodic convolution of regularly-spaced impulses with the DFT of my original sequence. Essentially, its like doing circular shifts on the input sequence, multiplying particular points by certain weights according to the square wave DFT, and summing them together. This can be represented as a linear system.

So I've generated the sparse matrix representing the convolution, and am now trying to solve Ax=b where A is the sparse matrix, x is the DFT of the original sequence, and b is the DFT of the corrupt sequence.

Unfortunately, the sparse matrix is ill conditioned. Using matlab, neither mldivide nor svds provided useful answers. The thing is, I know what the input sequence (and its DFT) look like. It's a DC component, two bin-centered sinusoids, and white noise. So I have an excellent initial guess, I just can't figure out a way to iterate on that. Also worth noting is that although the DFT is complex-valued, and these operations are all on complex numbers, I only care about the magnitudes.

If I solve mag(A)x=mag(b) I get a decent result, but its slightly too optimistic, so I'm hoping an iterative solution on the complex values will be more accurate.

My criteria for an acceptable solution is one that minimizes the mean square error between the magnitudes of Ax and b.

Any help would be greatly appreciated! I have little experience in these sorts of linear systems, but am doing my best to learn the appropriate methods.

This is the code I'm using to test my method: https://pastebin.com/F3bArgnR

• How ill-conditioned? (There's condest for checking.) Also, it sounds like you are better off solving a least squares problem instead; you could try lsqr or symmlq for that. (Both take initial guesses.) – Christian Clason Apr 30 '18 at 21:19
• condest = 9.4512e+16 That seems high, though honestly I have no idea what it actually means. lsqr seems to be promising, though. It's still a little off, but I suspect its because the system has a solution where a few elements are ~20 orders of magnitude larger than the rest, so I'm going to try splitting it up. It also only does one iteration. Is there a way to force more? Just crank down the tolerance? – Kayson Apr 30 '18 at 22:28
• Ok, so here's an interesting behavior - lsqr seems to work fine. And it gives me a valid solution to the system in that Ax=b to well within the accuracy I need. However, somehow the noise power in the calculated x is always lower than the noise power in the original x. And it is lower by almost exactly the factor 6/8. If I change the script such that I only have 5 of 8 samples, then that factor changes to 5/8. Same with 4 of 8 samples and 1/2... – Kayson Apr 30 '18 at 23:07
• That condition number shows that double precision will be devastated. Try 128 bit precision. – user14717 May 1 '18 at 10:05
• You could try regularizing for sparsity in the Fourier domain- this would get you back a signal with a small number of frequency components that matches the data reasonably well. – Brian Borchers May 1 '18 at 22:40