[SciPy-user] [Fwd: Re: Computing eigenvalues by trace minimization]

[dmitrey]

1) you don't need to use .T in dot, it's executing automatically:
f = lambda x: dot(x.T,dot(A,x)) # => dot(x,dot(A,x))
h = lambda x: dot(x.T,dot(B,x))-1.0 # => dot(x,dot(B,x))-1.0         
2) I have no symeig module, I will try to install it now
3) I didn't decide yet what to do if stop creteria funtol or xtol report 
true, but constraints are not satisfied yet (despite alg can 
successfully continue to decriase them, but I can't know will he do the 
trick or no and how many iterations it will require).
So now you should decide, what contol is ok for your problem, and then 
set appropriate contol, funtol, xtol
note that funtol and xtol are stop criteria, not desired tolerance of 
solution.
So, if you need contol=1e-6, as it is set by default, you need to reduce 
funtol and xtol (also, maybe, gradtol, if it will stop your problem 
before desired contol is achieved) from default values (1e-6) to 
something like 1e-8.
4) please wait some time before I will implement some changes about df, 
dc and dh. I will inform you about a hour later (as well as about my 
symeig module installation results) to update svn.
Use only func values for now.
5) Some weeks later (I hope) there will be excellent native OO QPQC 
solver (for positive-defined matrices only).

HTH, D.

Nils Wagner wrote:
>  Dmitrey,
>
> You asked me to move my inquiry to scipy-user.
> How can I improve my script wrt. the results ?
> How can I provide the gradients ?
> Is it possible to extend the code to rectangular matrices  instead of
> vectors ?
>
> python -i test_qpqc.py
> starting solver lincher (BSD  license)  with problem  unnamed
> itn 0: Fk= 0.285714285714 maxResidual= 0
> itn 10 : Fk= 0.0820810434319 maxResidual= 0.00260656014331
> N= 1.0 alpha= 0.978713763748
> itn 20 : Fk= 0.0810194693773 maxResidual= 1.71023808193e-05
> N= 1.0 alpha= 0.132923656348
> itn 22 : Fk= 0.0810154618808 maxResidual= 3.51005361554e-06
> N= 1.0 alpha= 0.354395906532
> solver lincher finished solving the problem unnamed
> istop:  4 (|| F[k] - F[k-1] || < funtol)
> Solver:   Time elapsed = 0.24   CPU Time Elapsed = 0.24
> NO FEASIBLE SOLUTION is obtained (max residual = 3.51005361554e-06)
>
> x_opt: [ 0.11973864  0.22980377  0.32143896  0.38732479  0.42178464 
> 0.4223829
>   0.38829845  0.32320301  0.2311393   0.12059048]
> f_opt: [ 0.08101546]
>
>
> Smallest eigenvalue by symeig 0.081014052771
>
> Nils
>
>   
> ------------------------------------------------------------------------
>
> from scikits.openopt import NLP
> from scipy import diag, ones, identity, rand, trace, dot, shape, zeros
> from scipy.linalg import solve, norm
> from symeig import symeig
>
> #f = lambda x: trace(dot(x.T,dot(A,x))) # Version for m > 1
> #h = lambda x: dot(x.T,dot(B,x))-identity(m)
> f = lambda x: dot(x.T,dot(A,x))
> h = lambda x: dot(x.T,dot(B,x))-1.0                 
>
> n = 10 # order of A and B
> m = 1  # subspace dimension # m > 1 doesn't work
>
> A = diag(2*ones(n))-diag(ones(n-1),-1)-diag(ones(n-1),1)
> B = identity(n)
>
> #x0 = rand(n,m)          # Initial subspace
> b = zeros(n)
> b[-1] = 1.
> x0 = solve(A,b)
> x0 = x0/norm(x0)
> p = NLP(f, x0, h=h)
>
> #Providing gradients is appriciated.
> #p.df = lambda x: 2*dot(A,x)
> #p.dh = lambda x: 2*dot(B,x)
>
> r = p.solve('lincher')
> print
> print 'x_opt:',r.xf
> print 'f_opt:',r.ff
> print
>
> w,v = symeig(A,B)
> print
> print 'Smallest eigenvalue by symeig',w[0]
>   
>
> ------------------------------------------------------------------------
>
> Subject:
> Re: Computing eigenvalues by trace minimization
> From:
> dmitrey <openopt at ukr.net>
> Date:
> Thu, 28 Jun 2007 17:48:21 +0300
> To:
> Nils Wagner <nwagner at iam.uni-stuttgart.de>
>
> To:
> Nils Wagner <nwagner at iam.uni-stuttgart.de>
>
>
> Hi Nils,
> your problem seems to be QPQC - quadratic problem with quadratic 
> constraints.
> Currently you can yuse only lincher to solve the problem.
> f = lambda x: trace (X^T A X)
>
> h = lambda x X^T B X - I_p
>
> p = NLP(f, x0, h=h)
> r = p.solve('lincher')
>
> Providing gradients is appriciated.
> p.df = ...
> p.dh = ...
>
> However, lincher can fail to solve the problem if something like 
> ill-conditioned matrix will be encountered.
> Our department's chiefman Petro I. Stetsyuk said me he can provide an 
> excellent QPQC algorithm for my Python QPQC solver (based on ralg), 
> but currently he is busy with other problems. Maybe, it will be ready 
> some weeks or months later.
> HTH, D.
>
>
>
> Nils Wagner wrote:
>> Hi Dmitrey,
>>
>> This inquiry is not related to my previous problem.
>> Can I use openopt to solve the quadratic minimization problem with 
>> openopt
>>
>> minimize trace (X^T A X)
>>
>> subjected to the constraints
>>
>> X^T B X = I_p,
>>
>> where I_p denotes the identity matrix of order p. A is symmetric
>> positive semidefinite
>> and B is symmetric positive definite. A, B \in \mathds{R}^{n \times n}.
>> A small example would be appreciated.
>>
>>
>>
>> Nils
>>  
>>
>> Reference: A. Sameh and J. Wisniewski
>> A trace minimization algorithm for the generalized eigenvalue problem
>> SIAM J. Numer. Anal. Vol. 19 (1982) pp. 1243-1259
>>
>>
>>
>>   
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>