[SciPy-user] Usage of fmin_ncg

[Nils Wagner]

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.


Nils

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)
    data.append(res)
    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 """
    return
dot(bsolve(v).T,v)*bsolve(p)+sigma_solve(p)+2*dot(outer(bsolve(v),bsolve(v)),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)*( )

random.seed(10)
v_0 = random.rand(n)                               # Initial vector

data=[]
v = optimize.fmin_ncg(F,v_0,Fp, fhess_p = Fhess_p)
semilogy(arange(0,len(data)),data)
print  'Newton CG',-1./(2*sqrt(-F(v)))+sigma
print R(bsolve(v))
print