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@ukr.net> Date: Thu, 28 Jun 2007 17:48:21 +0300 To: Nils Wagner <nwagner@iam.uni-stuttgart.de>

To: Nils Wagner <nwagner@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@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user