One of my favorite formulas is the factorization of :
This is a fairly ubiquitous formula that most people would have seen in high school, perhaps in various guises. A common special case is which is itself a special case of the “sum of squares” formula
The irreducible factors of are the cyclotomic polynomials which have many interesting links to number theory.
In other areas, substituting something nice for leads to interesting formulas/results. When is some prime power (usually written ), the quantity
is the -analog of The reason for the name is that as we let we get
However, I won’t be writing about that today. Instead, this post is about a curious connection between that simple formula above and deconvolution.
Geometric series and the Neumann series
For this purpose, let’s start by letting be small (e.g. ), then so
Letting , we obtain the geometric series which is the Taylor series expansion for
Going one step further, replacing with a matrix we get the Neumann series
If is suitably small in some sense (so that for large enough ), then this infinite series converges to the inverse of where is the identity matrix. By truncating the series, we get
which approximates for suitably large If, further, is nilpotent (i.e. for some ) then is the exact inverse of
The Neumann series thus provides an iterative way of solving linear equations if the system can be written as for some suitably small If we have the equation
and we’re trying to solve for given , define a recursive relation:
Observe that this is equivalent to the recursion
which makes it obvious that
When is large enough, will be small, and will approximate the solution
What has all this got to do with deconvolution?
Now suppose I want to record a signal (e.g. photo) with some device (e.g. camera). For simplicity, let’s say my signal is just an array of values:
Unfortunately, during my measurement process, my signal gets “blurred”, causing me to record this instead:
In technical terms, my signal is a convolution of the original signal with some point-spread function. To recover the true green signal from the corrupted blue one, I’ll need to carry out deconvolution.
But it turns out that in this (grossly oversimplified) case, we already know how to deconvolve this signal!
Let’s say I have some additional information about the true signal: for example, I may expect the original signal to have only one peak, and it seems like the corrupted signal has 2 peaks. I may thus guess that my recorded signal is the sum of the original signal with a translated (and possibly scaled) version of the original signal such as this:
As you can see, it is indeed true that the blue signal = green + red signals, and the translation is 2 units to the right. Of course, if I know the red signal, then I can subtract it from the blue signal to get the green signal. But what if all I have is the blue signal?
If we treat the (unknown) green signal as a vector , then translation to the right by 2 units is given by this matrix
This is a nilpotent matrix, because , as can be easily verified by repeated matrix multiplication, or by the observation that if we translate 2 units to the right 4 times, we’ll just end up with all zeros.
So we guess that our blue signal is given by , where is the identity matrix and is some scaling constant that I have to guess. Letting , we can recover our original signal via the iterative procedure above:
At each iteration, we’re adding to . Here’s what each looks like; they are essentially scaled, translated copies of the blue signal:
It might seem strange that to recover the original signal, we take even more corrupted copies of the corrupted signal. But this is precisely what’s going on in the formula
The formula says that to invert , we just sum many copies (infinitely many, too!) of powers of !
Back to reality
In this case, we obtained perfect reconstruction of our original green signal. This was only possible because I used the correct value of . Here’s what would’ve happened if I didn’t know the true value of :
It’s pretty close, but not a perfect reconstruction of the original signal.
“Naïve deconvolution” as I described above looks good on paper, but will seldom work with real data. Recall that to successfully get the green signal, I had to
- Guess that the corrupted data was due to a simple translation
- Guess that the translation by was 2 units to the right (by looking at the peaks)
- Guess the correct scaling factor
Essentially, I had to successfully guess that the point-spread function looked like this:
This is about as simple as point-spread functions go. If signals were photos, the blue signal might look like this:
If you look closely, you can see that the stars are “doubled”. This was caused by me accidentally knocking into my camera tripod while the photo was being taken.
The double stars give a pretty good gauge of what the corrupting transformation was: a very slight rotation, centered somewhere above the top-left corner of the photo. Even knowing this, I had to carry out many rounds of trial and error to find the rotation parameters, as well as the intensity scaling factor (the double stars are about equally bright, so ), and even then, as the residual streaks at the top right corner of the next image indicate, my parameters are still off.
In the deconvoluted image, the stars have been contracted to points (not quite perfectly). But there are lots of “wrinkles” in the image. This is because instead of obtaining
what I got was a noisy version . At each iteration of the deconvolution procedure, the noise gets multiplied and propagated by !
There are also vertical bars at the bottom left of the image. These arise because of the artificial truncation that occurs at the boundary of the image. To remove these bars, one would have to guess what the image looks like outside the photo.
Even in the absence of noise or truncation, it is often the case that the point-spread function results in an that is not “small enough”. Recall that for our method to work, we needed as . However, many point-spread functions (such as Gaussian blur) result in matrices that blow up as . Using this iterative method will only make the image worse.
Difficulties like these are why people can spend their lives studying various methods of deconvolution. But if you can deconvolve well, you’re practically guaranteed a job in image processing, finance, microscopy, seismology and basically any field that requires signal processing!
Not great, but not too shabby for a humble high school formula!