[SciPy-dev] Trouble with fmin_ncg

[dmitrey]

Nils, I'm not sure that developers related to scipy.sparse read the 
messages with "Trouble with fmin_ncg" topic.
You better to rename the one and/or call them directly.
HTH, D.

Nils Wagner wrote:
> Hi,
>
> I guess I traced the problem with ncg down to the Hessian.
> The results between the sparse and the dense version differ.
>
> Nils
>  
>   
> ------------------------------------------------------------------------
>
> from numpy import dot,  outer, random, argsort, arange, array
> from scipy import io, linalg, optimize
> from scipy.sparse import speye
>  
> def Rd(v):
>     """ Rayleigh quotient with dense arrays """
>     rq = dot(v,dot(A1,v))/dot(v,dot(B1,v))
>     return rq                                 
>
> def Rs(v):
>     """ Rayleigh quotient with sparse matrices """
>     rq = dot(v,A*v)/dot(v,B*v)
>     return rq
>
> def Rppd(v):
>     """ Hessian with dense arrays  """
>     result = 2*(A1-Rd(v)*B1-outer(dot(B1,v),Rp(v))-outer(Rp(v),dot(B1,v)))/dot(v.T,dot(B1,v))
>     return result
>
> def Rpps(v):
>     """ Hessian with sparse matrices"""
>     result = 2*(A-Rs(v)*B-outer(B*v,Rps(v))-outer(Rps(v),B*v))/dot(v,B*v)
>     return result
>
> def Rp(v):
>     """ Gradient with dense arrays"""
>     result = 2*(dot(A1,v)-Rd(v)*dot(B1,v))/dot(v.T,dot(B1,v))
>     return result
>
> def Rps(v):
>     """ Gradient with sparse matrices """
>     result = 2*(A*v-Rs(v)*B*v)/dot(v.T,B*v)
>     return result
>                                         
> A = io.mmread('nos4.mtx') 
> n = A.shape[0]
> A1 = A.todense()
>
> B = speye(n,n)
>
> B1 = B.todense()
>
> A1 = array(A1)
> B1 = array(B1)
>
> for i in arange(0,2):
>     
>     v_0=random.rand(n)
>     print
>     print
>     print
>     print Rd(v_0)-Rs(v_0)
>     print linalg.norm(Rp(v_0)-Rps(v_0))
>     print linalg.norm(array(Rpps(v_0)) - Rppd(v_0))
>     
>
>   
> ------------------------------------------------------------------------
>
> _______________________________________________
> Scipy-dev mailing list
> Scipy-dev at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-dev
>