[SciPy-User] Request help with fsolve outputs
Moore, Eric (NIH/NIDDK) [F]
eric.moore2 at nih.gov
Fri Oct 26 16:53:47 EDT 2012
> -----Original Message-----
> From: The Helmbolds [mailto:helmrp at yahoo.com]
> Sent: Friday, October 26, 2012 11:52 AM
> To: User SciPy
> Subject: [SciPy-User] Request help with fsolve outputs
>
> 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[0], x[1], x[2]).
>
> 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?
>
>
> Bob H
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
I haven't spent any time playing with your example, but I have looked at a simpler example. It seems that the approximated jacobian can be quite far off and fsolve can still return correct zeros. How good the approximation ends up being depends on your choice of initial conditions. I've attached a trivial example that shows this.
There is also the epsfcn parameter, but based on the minpack documentation, the default should be okay for both our examples since we have full machine precision available in our functions to minimize.
It would make things much easier for other people if you could post a python file where you have defined your function and are calling fsolve. I believe Ralph asked if you could do this the last time you posted this question. I'm sorry that no one has come in with a simple answer, but this would still be a good step to make considering your question a little easier for other people.
Good luck,
Eric
-------------- next part --------------
A non-text attachment was scrubbed...
Name: fsolve_test.py
Type: application/octet-stream
Size: 583 bytes
Desc: fsolve_test.py
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20121026/ffda53ba/attachment.obj>
More information about the SciPy-User
mailing list