[SciPy-User] Request help with fsolve outputs
Fri Oct 26 10:51:36 CDT 2012
Please help me out here. I’m trying to rewrite the docstring for the `fsolve.py` routine
located on my machine in: C:/users/owner/scipy/scipy/optimize/minpack.py
The specific issue I’m having difficulty with is understanding the outputs described in fsolve’s docstring as:
'fjac': the orthogonal matrix, q, produced by the QR factorization of the final approximate Jacobian matrix, stored column wise.
'r': upper triangular matrix produced by QR factorization of same matrix.
These are described in SciPy’s minpack/hybrd.f file as:
‘fjac’ is an output n by n array which contains the orthogonal matrix q produced by the qr factorization of the final approximate jacobian.
‘r’ is an output array of length lr which contains the upper triangular matrix produced by the qr factorization of the final approximate jacobian, stored rowwise.
For ease in writing, in what follows let’s use the symbols ‘Jend’ for the final approximate Jacobian matrix, and use ‘Q’ and ‘R’ for its QR decomposition matrices.
Now consider the problem of finding the solution to the following three nonlinear
equations (which we will refer to as 'E'), in three unknowns (u, v, w):
2 * a * u + b * v + d - w * v = 0
b * u + 2 * c * v + e - w * u = 0
-u * v + f = 0
where (a, b, c, d, e, f ) = (2, 3, 7, 8, 9, 2). For inputs to fsolve, we identify
(u, v, w) = (x, x, x).
Now fsolve gives the solution array:
[uend vend wend] = [ 1.79838825 1.11210691 16.66195357].
With these values, the above three equations E are satisfied to an accuracy of about 9 significant figures.
The Jacobian matrix for the three LHS functions in E is:
J = np.matrix([[2*a, b-w, -v], [b-w, 2*c, -u], [-v, -u, 0.]])
Note that it’s symmetrical, and if we compute its value using the above fsolve’s ‘end’ solution values we get:
Jend = [[ 4. 19.66195357 1.11210691],
[ 19.66195357 14. 1.79838825],
[ 1.11210691 1.79838825 0. ]]
Using SciPy’s linalg package, this Jend has the QR decomposition:
Qend = [[-0.28013447 -0.91516674 -0.28981807]
[ 0.95679602 -0.24168763 -0.16164302]
[ 0.07788487 -0.32257856 0.94333293]]
Rend = [[-14.278857 17.08226116 -1.40915124]
[ -0. 9.69946027 1.45241144]
[ -0. 0. 0.61300558]]
and Qend * Rend = Jend to within about 15 significant figures.
However, fsolve gives the QR decomposition:
qretm = [[-0.64093238 0.75748326 0.1241966 ]
[-0.62403598 -0.60841098 0.4903215 ]
[-0.44697291 -0.23675978 -0.8626471 ]]
rret = [ -7.77806716 30.02199802 -0.819055 -10.74878184 2.00090268 1.02706198]
and converting rret to an upper triangular NumPy matrix gives:
rretm = [[ -7.77806716 30.02199802 -0.819055 ]
[ 0. -10.74878184 2.00090268]
[ 0. 0. 1.02706198]]
Now qret and rretm bear no obvious relation to Qend and Rend.
Although qretm is orthogonal to about 16 significant figures, we find the product:
qretm * rretm = [[ 4.98521509 -27.38409295 2.16816676]
[ 4.85379376 -12.19513008 -0.2026608 ]
[ 3.47658529 -10.87414051 -0.99362993]]
which bears no obvious relationship to Jend.
The hybrd.f routine in minpack refers to a permutation matrix, p, such that we should have in our notation:
p*Jend = qretm*rretm,
but fsolve apparently does not return the matrix p, and I don’t see any permutation of Jend that would equal qretm*rretm.
The hybrd.f routine does refer to some "scaling" that is going on, but my Fortran is about 40 years too stale for me to interpret it.
If we reinterpret rret as meaning the matrix:
rretaltm = [[ -7.77806716 30.02199802 -10.74878184]
[ 0. -0.819055 2.00090268]
[ 0. 0. 1.02706198]]
then we get the product:
qretm * rretaltm = [[ 4.98521509 -19.86249109 8.53245022]
[ 4.85379376 -18.2364849 5.99384603]
[ 3.47658529 -13.22510045 3.44468895]]
which again bears no obvious relationship to Jend. Using the transpose of qretm in the above product is no help.
So please help me out here. What are the fjac and r values that fsolve returns?
How are they related to the above Qend, Rend, and Jend?
How is the user supposed to use them?
More information about the SciPy-User