Corrected mistakes in functioning of SVD
Corrected mistakes in functioning of SVD

Here is how SVD functions as of this commit.

Let's say we have 2 parameters and 20 data points. The Jacobian is
[20, 2] and the residual vector is [20, 1].

We then used reduced SVD (not full SVD, see the NumPy documentation
for numpy.linalg.svd) on the Jacobian, J = U . s . VT.

Here, U is [20, 2], s is a [2, 1] vector and VT is the transpose of
V, both of which are [2, 2]. In code, this is as so.

mu, vs, mvt = np.linalg.svd(jacob)
ms = np.diag(vs)

We can form the diagonal matrix, ms or S, from vs or s as shown
above. We can also confirm that this method works as intended by
doing the following.

reform = mu.dot(np.diag(vs).dot(mvt))
np.allclose(jacob, reform)

This should return True, indicating that jacob and reform are the
same. This check is now built into the function return_svd, but is
not used by default in order to save time.

We can form the inverted matrix, S', as so.

ms = np.diag(vs)
msi = np.linalg.inv(ms)

In this inverted matrix, msi or S', smaller values come 1st (diagonal
elements increase from index [0, 0] to [n, n] where n is the length
of vector vs).

~~~~~~~~~~

SVD without any modifications is performed to generate trial parameter
changes as such.

X' = VT . S' . UT . r
X' = [2, 2] . [2, 2] . [2, 20] . [20, 1]
X' = [2, 2] . [2, 2] . [2, 1]
X' = [2, 2] . [2, 1]
X' = [2, 1]

Here, X' is the parameter changes and r is the residual vector.

In code, here's what it looks like.

changes = mvt.T.dot(msi.dot(mu.T.dot(resid)))

~~~~~~~~~~

We then modify S' in various way to acquire other trial parameter
changes.

For SVD without thresholds, we zero diagonal elements of S'/msi to
form S", starting with the largest values. A simplified version of
the code looks something like this.

for i in xrange(0, len(vs) - 1):
    msi[-(i + 1), -(i + 1)] = 0.
    changes = mvt.T.dot(msi.dot(mu.T.dot(resid)))

Here is an example of how S"/msi evolves over time for a system
involving 4 parameters.

------------------------------ ZEROED 0 ELEMENTS ------------------------------
>>> msi:
[[  6.20209407e-01   0.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  0.00000000e+00   7.95221576e-01   0.00000000e+00   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   1.82438431e+00   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00   1.24827979e+16]]
>>> changes:
[ -1.49218287e+16   1.63099821e+16   2.54761516e+00  -1.01361975e+16]
------------------------------ ZEROED 1 ELEMENTS ------------------------------
>>> msi:
[[ 0.62020941  0.          0.          0.        ]
 [ 0.          0.79522158  0.          0.        ]
 [ 0.          0.          1.82438431  0.        ]
 [ 0.          0.          0.          0.        ]]
>>> changes:
[ 2.98237556  2.05417715  1.1976335  -1.08511153]
------------------------------ ZEROED 2 ELEMENTS ------------------------------
>>> msi:
[[ 0.62020941  0.          0.          0.        ]
 [ 0.          0.79522158  0.          0.        ]
 [ 0.          0.          0.          0.        ]
 [ 0.          0.          0.          0.        ]]
>>> changes:
[ 2.86364141  1.94166881  1.55287877 -1.09135434]
------------------------------ ZEROED 3 ELEMENTS ------------------------------
>>> msi:
[[ 0.62020941  0.          0.          0.        ]
 [ 0.          0.          0.          0.        ]
 [ 0.          0.          0.          0.        ]
 [ 0.          0.          0.          0.        ]]
>>> changes:
[ 0.52200671 -0.37786601  0.03060935 -1.37648095]

~~~~~~~~~~

For SVD with thresholds, the user supplied thresholds are sorted in reverse.

self.svd_factors = [0.001, 0.01, 0.1, 1.]
self.svd_factors.sort(reverse=True)
self.svd_factors
> [1., 0.1, 0.01, 0.001]

Using this method, S"/msi evolves as like so.

---------------------------------- NO FACTORS ---------------------------------
>>> msi:
[[  6.20209407e-01   0.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  0.00000000e+00   7.95221576e-01   0.00000000e+00   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   1.82438431e+00   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00   1.24827979e+16]]
>>> changes:
[ -1.49218287e+16   1.63099821e+16   2.54761516e+00  -1.01361975e+16]
--------------------------------- FACTOR: 1.0 ---------------------------------
>>> msi:
[[ 0.62020941  0.          0.          0.        ]
 [ 0.          0.79522158  0.          0.        ]
 [ 0.          0.          0.          0.        ]
 [ 0.          0.          0.          0.        ]]
>>> old_msi:
[[  6.20209407e-01   0.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  0.00000000e+00   7.95221576e-01   0.00000000e+00   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   1.82438431e+00   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00   1.24827979e+16]]
>>> changes:
[ 2.86364141  1.94166881  1.55287877 -1.09135434]
--------------------------------- FACTOR: 0.1 ---------------------------------
>>> msi:
[[ 0.  0.  0.  0.]
 [ 0.  0.  0.  0.]
 [ 0.  0.  0.  0.]
 [ 0.  0.  0.  0.]]
>>> old_msi:
[[ 0.62020941  0.          0.          0.        ]
 [ 0.          0.79522158  0.          0.        ]
 [ 0.          0.          0.          0.        ]
 [ 0.          0.          0.          0.        ]]
  -- Vector is all zeros. Breaking.

~~~~~~~~~~

Also, addressed issue #14 of q2mm/q2mm (not this repository). This
issue stated:

"""
This concerns using SVD without thresholds.

To recap, J = UsV.

This method sequentially zeros elements from s. However, I think it
goes one step too far and zeros every last element of s, so of course
the parameter changes end up being zero across the board.
"""

Now the code stops before having matrix S being entirely zero.

ericchansen/q2mm,arosale4/q2mm,Q2MM/q2mm,Q2MM/q2mm,arosale4/q2mm,ericchansen/q2mm