note:
examples of the matrices can be found in
demo5.html
demo5b.html
I perform a blurring of a vector b by averaging neighbouring elements.
For this, I use the matrix H:
|1/h 1/h 1/h 1/h 1/h 0 0 0 0 0 |
| 0 1/h 1/h 1/h 1/h 1/h 0 0 0 0 |
| 0 0 1/h 1/h 1/h 1/h 1/h 0 0 0 |
| 0 0 0 1/h 1/h 1/h 1/h 1/h 0 0 |
| 0 0 0 0 1/h 1/h 1/h 1/h 1/h 0 |
| 0 0 0 0 0 1/h 1/h 1/h 1/h 1/h |
I call H the "blurring matrix".
So
H b = x, and I want to find b from x and H.
The number of non-zero elements on a row is h, which I call the "spread".
The number of columns equals the length of b, which I will call N.
the number of rows = N-h.
I will write g for my estimate of b, so
H g = x
Note that H is not square, and there may be several (actually, infinitely
many) solutions for g (because H has more columns than rows).
I add a few rows of a unit matrix, to obtain an invertible matrix:
A =
| 1 0 0 0 0 0 0 0 0 0 |
| 0 1 0 0 0 0 0 0 0 0 |
| 0 0 1 0 0 0 0 0 0 0 |
| 0 0 0 1 0 0 0 0 0 0 |
|1/h 1/h 1/h 1/h 1/h 0 0 0 0 0 |
| 0 1/h 1/h 1/h 1/h 1/h 0 0 0 0 |
| 0 0 1/h 1/h 1/h 1/h 1/h 0 0 0 |
| 0 0 0 1/h 1/h 1/h 1/h 1/h 0 0 |
| 0 0 0 0 1/h 1/h 1/h 1/h 1/h 0 |
| 0 0 0 0 0 1/h 1/h 1/h 1/h 1/h |
and I extend the vector x in a clever way, and then use the inverse of A to find g.
So I define y = [ c0, c1, c2, ..., c_h ,x0, x1, .... , x_n].
We then have
| 1 0 0 0 0 0 0 0 0 0 | |g0| |c0|
| 0 1 0 0 0 0 0 0 0 0 | |g1| |c1|
| 0 0 1 0 0 0 0 0 0 0 | |g2| |c2|
| 0 0 0 1 0 0 0 0 0 0 | |g3| |c3|
|1/h 1/h 1/h 1/h 1/h 0 0 0 0 0 | |g4| = |x0|
| 0 1/h 1/h 1/h 1/h 1/h 0 0 0 0 | |g5| |x1|
| 0 0 1/h 1/h 1/h 1/h 1/h 0 0 0 | |g6| |x2|
| 0 0 0 1/h 1/h 1/h 1/h 1/h 0 0 | |g7| |x3|
| 0 0 0 0 1/h 1/h 1/h 1/h 1/h 0 | |g8| |x4|
| 0 0 0 0 0 1/h 1/h 1/h 1/h 1/h | |g9| |x5|
The effect of the extension of x
The inverse looks like:
inv(A) =
| 1 0 0 0 0 0 0 0 0 0 0 0 0 |
| 0 1 0 0 0 0 0 0 0 0 0 0 0 |
| 0 0 1 0 0 0 0 0 0 0 0 0 0 |
| 0 0 0 1 0 0 0 0 0 0 0 0 0 |
|-1 -1 -1 -1 h 0 0 0 0 0 0 0 0 |
| 1 0 0 0 -h h 0 0 0 0 0 0 0 |
| 0 1 0 0 0 -h h 0 0 0 0 0 0 |
| 0 0 1 0 0 0 -h h 0 0 0 0 0 |
| 0 0 0 1 0 0 0 -h h 0 0 0 0 |
|-1 -1 -1 -1 h 0 0 0 -h h 0 0 0 |
| 1 0 0 0 -h h 0 0 0 -h h 0 etc.... |
So, the inverse is a triangular matrix, in which a small unit matrix is repeated in the first columns.
So, when the multiplication
g= inv(A) y
is performed, each offset c0,...,c_h will have effect on several, regularly spaced, elements of g.
If I choose one offset very high, it will also cause a repeated error. I call such an error a "ripple" because it repeats with a fixed period.
But if I make all offsets 0, it is very likely that repeated occurences in the same columns of "-h h" will also lead to a ripple (repeated error).
The rows with "-1 -1 -1 -1 ..." mean that the total value of the offsets will cause a ripple if their sum is too big.
So I want to choose c0,...,c_h such that, on average, they don't have effect.
For this reason, I want to extract the average value of each ripple.
So I take "periodic subsets" of g (where the period is h, of course):
C g = C inv(A) y
where
C =
|1/m 0 0 0 0 1/m 0 0 0 0 1/m 0 0 0 0 1/m 0 |
| 0 1/m 0 0 0 0 1/m 0 0 0 0 1/m 0 0 0 0 1/m|
| 0 0 1/n 0 0 0 0 1/n 0 0 0 0 1/n 0 0 0 0 |
| 0 0 0 1/n 0 0 0 0 1/n 0 0 0 0 1/n 0 0 0 |
| 0 0 0 0 1/n 0 0 0 0 1/n 0 0 0 0 1/n 0 0 |
In this example, m=4 and n=3.
The number n or m equals the number of non-zero elements in the row. In the
first rows, there are 4 non-zero elements, so m=4. in the last rows, there
are 3 (in this example).
I call C the "mask", and C inv(A) the "calibration Matrix"
Now, C inv(A) looks someting like
|1 0 0 0 -3.33 3.33 0 0 0 -1.67 1.67|
|0 1 0 0 0 -2.5 2.5 0 0 0 0 |
|0 0 1 0 0 0 -2.5 2.5 0 0 0 |
|0 0 0 1 0 0 0 -2.5 2.5 0 0 |
|-1 -1 -1 -1 5 0 0 0 -2.5 2.5 0 |
In this example, I have 4 offsets (c0,...,c3) to adjust the result C inv(A) y
which is a vector of 5 elements!
(if you look at the definition of matrix A, you will notice that I added 4 offsets, when the blurring involves 5 elements).
To start, I separate the effects of (c0,..,c_h) and (x0,..., xN).
To this effect, I separate the "calibration matrix" into the "offsets Matrix" F and the "ripple Matrix" R:
| F | | R | =
|1 0 0 0 | |-3.33 3.33 0 0 0 -1.67 1.67|
|0 1 0 0 | |0 -2.5 2.5 0 0 0 0 |
|0 0 1 0 | |0 0 -2.5 2.5 0 0 0 |
|0 0 0 1 | |0 0 0 -2.5 2.5 0 0 |
|-1 -1 -1 -1 | |5 0 0 0 -2.5 2.5 0 |
It would be nice to obtain:
F c + R x = 0
by a proper choice of c. Of course, this might be impossible, because F has more rows (equations to satisfy) then columns (variables).
My intuitive thought was then, to calculate the least-squares approximation for c, and there is a good reason why this is a valid approach.
Actually, an offset for all ripples is no problem, so it is sufficient to require
F c + R x = r [1, 1, 1, ... , 1] (with r some real number), because:
If (R x) looks like this (in a graph)
value -->
___________________
index | *
| | *
V | *
| *
| *
(-F c) may look like this:
value -->
___________________
index | %
| | %
V | %
| %
| %
because the difference (R x - F c)[i] is the same for each element.(this is clear if I plot both in 1 graph).
value -->
___________________
index | * %
| | * %
V | * %
| * %
| * %
And for the matrix F, this is obtained by the least squares approximation!
(because an increase of a offset c_i will increase (F c)[i] with the same amount, and decrease the last element (F c)[h]).
The squares will be minimal if
for each index i: (F c - R x)[i] = (F c - R x)[h].
Therefor, I calculate the least-squares approximation of c. (I will not explain that procedure here, it is public knowledge).
Are the offsets part of the solution g?
Yes.
Every element of x contains information about an element of b, and its neighbours in b.
(this is a consequence of the matrix H).
So it is likely that when I unblur, I obtain information about elements that are outside x.
This is valid at the start and at the end of x.
The offsets are calculated with the objective to minimize ripples.
A ripple occurs if the difference between two elements of x is not adequately compensated (by the values of other elements and ultimately, by the offsets).
(a look at inv(A) may help to understand this).
Actually, the removal of ripples is achieved by the fact that the values of the first elements of g are fitting with the values of x.
In other words, the elimination of ripples forces the offsets to resemble the original values.
Note that the solution is an approximation, or even: a compromise.
Elements of x that have a "median" index are not treated with more or less preference than the first or last elements.
So, for the solution g, the values for the last indices should be as reliable as for the first indices.
The offsets are not calculated by inversion of A. Do they approach the original values of b as well as the elements at the end of g do?
Numerical verification suggests that the offsets are calculated as precisely as other elements of g. But the reason is:
ALL elements of g are calculated to satisfy
H g = x
under the condition that:
"the average of each periodic subset has the same value". (otherwise, there is a "repeated deviation" or ripple).
It will not make a difference whether I choose for the offsets elements at the start, at the end, or somewhere in between (as long as they are not on the same "ripple" of course).
Numerical verification
Execute "demo5.html".
Note that an internet connection is required for the best rendering of the matrices.
I use the following variable names for variables:
b:B
g:foundB
y: foundX
x: xTrunc
F: offsetsMatrix
R: rippleMatrix
Rx:ripple
-c: offsets (so it must be multiplied by -1 to obtain c)
calibration Matrix: calMat
And I perform some checks:
-Fc: approxRipple
C inv(A) y = calMat.foundX = calX
It turns out that calX = r [1, 1, .... , 1] with r a real number. So each
ripple has the same average.
The quality of the result (foundB) depends on:
- the original B
- the "spread" of the blurring (variable h)
- the length of B (the best is: long, compared to h)
Observations
1.
The reconstructed vector, foundB, is often better than the blurred vector, xTrunc.
Although it seems that, if the h is big compared to N, the suppression of the ripples is too strong: the offsets correspond well with the orignal, but there happens to be a ripple at other places.
2.
If the original (B) is symmetric (i.e. if B[i]=B[dim-i-1], foundB is symmetric too.
===> So the calculated offsets are part of the solution, in the same way as the last elements of the solution. <===
3. The desire to avoid ripples can lead to an error if elements with higher indices have, systematically, a higher value.
e.g. if B=[1 2 3 4 5 6 7 8 ..... N], the solution will look like
[ 0, 2.5, 2.5, 2.5, 2.5, 5, 7.5, 7.5, 7.5, 7.5, 10,..]
4. demo5b.html demonstrates that the results are the same if I place the offsets at a different position.
For example, by taking
A =
| 0 0 0 1 0 0 0 0 0 0 |
| 0 0 0 0 1 0 0 0 0 0 |
| 0 0 0 0 0 1 0 0 0 0 |
| 0 0 0 0 0 0 1 0 0 0 |
|1/h 1/h 1/h 1/h 1/h 0 0 0 0 0 |
| 0 1/h 1/h 1/h 1/h 1/h 0 0 0 0 |
| 0 0 1/h 1/h 1/h 1/h 1/h 0 0 0 |
| 0 0 0 1/h 1/h 1/h 1/h 1/h 0 0 |
| 0 0 0 0 1/h 1/h 1/h 1/h 1/h 0 |
| 0 0 0 0 0 1/h 1/h 1/h 1/h 1/h |