[SciPy-user] Usage of fmin_ncg

Nils Wagner nwagner@iam.uni-stuttgart...
Thu Mar 8 10:06:38 CST 2007

Hi all,

I have a question wrt optimize.fmin_ncg(F,v_0,Fp, fhess_p = Fhess_p).
especially the last argument is of interest. I am not sure if I applied
it correctly.

The help function yields

      fhess_p -- a function to compute the Hessian of f times an
                 arbitrary vector: fhess_p (x, p, *args)

Is p the arbitrary vector ?
Any hint will be appreciated.


P.S.: The Hessian is not directly available.

from scipy import *
from pylab import plot, show, semilogy, xlabel, ylabel, legend, savefig,
title, figure, ylim
def R(v):
    """ Rayleigh quotient """
    return dot(v.T,A*v)/dot(v.T,B*v)

def F(v):
    """ Objective function """
    Bv = bsolve(v)
    res = linalg.norm(A*Bv-R(Bv)*B*Bv)/linalg.norm(B*Bv)
    h =    0.25*dot(bsolve(v).T,v)**2+0.5*dot(sigma_solve(v).T,v)
    return h

def Fp(v):
    """ Gradient of the objective function """
    return dot(bsolve(v).T,v)*bsolve(v)+sigma_solve(v)

def Fhess_p(v,p):
    """ Hessian applied to an arbitrary vector p """

# Test matrices from MatrixMarket
A = io.mmread('bcsstk02.mtx.gz')
B = io.mmread('bcsstm02.mtx.gz')
n = A.shape[0]

sigma = 47.5                                       # Preassigned shift
bsolve = linsolve.splu(B).solve                    # inv(B)*( )
sigma_solve = linsolve.splu(A - sigma*B).solve     # inv(A-sigma*B)*( )

v_0 = random.rand(n)                               # Initial vector

v = optimize.fmin_ncg(F,v_0,Fp, fhess_p = Fhess_p)
print  'Newton CG',-1./(2*sqrt(-F(v)))+sigma
print R(bsolve(v))

More information about the SciPy-user mailing list