[SciPy-dev] Trouble with optimize.fmin_ncg

[dmitrey]

you have line
A = io.mmread('nos4.mtx') # clustered eigenvalues
but you didn't provide nos4.mtx file
please send it to me.
Regards, D.


Nils Wagner wrote:
> Hi all,
>
> I have some trouble with fmin_ncg. IIRC the attached code worked for me
> before (line 44-48) but it hangs with recent svn.
>
> Any pointer would be appreciated.
>  
> Nils
>
>   
> ------------------------------------------------------------------------
>
> from scipy import *
> from scipy.sparse import speye
> from pylab import plot, show, semilogy, xlabel, ylabel, figure, savefig, legend
>
> def R(v):
>     rq = dot(v.T,A*v)/dot(v.T,B*v)
>     res = (A*v-rq*B*v)/linalg.norm(B*v)
>     data.append(linalg.norm(res))
>     return rq                                 
>
> def Rp(v):
>     """ Gradient """
>     return 2*(A*v-R(v)*B*v)/dot(v.T,B*v)                                                              
>
> def Rpp(v):
>     """ Hessian """
>     return 2*(A-R(v)*B-outer(B*v,Rp(v))-outer(Rp(v),B*v))/dot(v.T,B*v)
>
>
> A = io.mmread('nos4.mtx') # clustered eigenvalues
> #B = io.mmread('bcsstm02.mtx.gz') 
> #A = io.mmread('bcsstk06.mtx.gz') # clustered eigenvalues
> #B = io.mmread('bcsstm06.mtx.gz') 
> n = A.shape[0]
> B = speye(n,n)
> random.seed(1)
> v_0=random.rand(n)
>
> if n < 1000:
>
>    w,vr = linalg.eig(A.todense(),B.todense())
>    ind = argsort(w.real)
>    print 'Least eigenvalue',w[ind[0]]    
>
>
> data=[]
> v,fopt, gopt, Hopt, func_calls, grad_calls, warnflag,allvecs = optimize.fmin_bfgs(R,v_0,fprime=Rp,full_output=1,retall=1)
> if warnflag == 0:
>    semilogy(arange(0,len(data)),data)
>    print 'Rayleigh quotient BFGS',R(v)
> #
> # The program hangs if fmin_ncg is active
> #
> #data=[]
> #v,fopt, fcalls, gcalls, hcalls, warnflag,allvecs = optimize.fmin_ncg(R,v_0,fprime=Rp,fhess=Rpp,full_output=1,retall=1)
> #if warnflag==0:
> #   semilogy(arange(0,len(data)),data)
> #   print 'Rayleigh quotient NCG',R(v)
>
> data=[]
> v,fopt, func_calls, grad_calls, warnflag,allvecs = optimize.fmin_cg(R,v_0,fprime=Rp,full_output=1,retall=1)
> semilogy(arange(0,len(data)),data)
> print 'Rayleigh quotient CG ',R(v),fopt
>
> data=[]
> rc,nfeval,v = optimize.fmin_tnc(R,list(v_0),fprime=Rp)
> #print optimize.tnc.RCSTRINGS[rc]
> v = array(v)
> semilogy(arange(0,len(data)),data)
> print 'Rayleigh quotient',R(v)
>
> data=[]
> v,f,d = optimize.fmin_l_bfgs_b(R,v_0,fprime=Rp)
> semilogy(arange(0,len(data)),data)
> print 'Rayleigh quotient l_bfgs_b',R(v),f                  
>
> xlabel('Function evaluation')
> ylabel(r'Residual $\|Ax-\lambda\,Bx\|/\|Bx\|$')
> #legend(('fmin\_bfgs','fmin\_ncg','fmin\_tnc','fmin\_l\_bfgs\_b'),shadow=True)
> legend(('fmin\_bfgs','fmin\_cg','fmin\_tnc','fmin\_l\_bfgs\_b'),shadow=True)
> #legend(('fmin\_bfgs','fmin\_ncg','fmin\_tnc','fmin\_l\_bfgs\_b'),shadow=True)
> #savefig('leastR.png')
>
> figure(2)
>
> semilogy([0],[R(v)],'ro')
> semilogy(arange(0,n),w[ind].real)
> xlabel(r'$i$')
> ylabel(r'$\lambda_i$')
> legend((r'least eigenvalue $\lambda_1$','Spectrum'),shadow=True,loc=4)
>
> show()
>   
> ------------------------------------------------------------------------
>
> _______________________________________________
> Scipy-dev mailing list
> Scipy-dev at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-dev
>