What type of optimization problem is this?

Specific Domains Optimization (Mathematical)
1

Hi all. I have a constrained optimization problem which looks like it might belong to a known class.

The problem is to maximize the quotient \dfrac{x^T A x}{x^T B x} w.r.t real vector x, subject to the constraint x^T x = 1. Matrices A, B are real n \times n and B is also symmetric positive definite. In general A is neither symmetric nor positive definite.

Do specialized algorithms exist for this type of problem? I suspect that I may need to use generic non-linear optimization. In either case, I’m interested to know whether solvers are available in Julia which might be effective for this problem. Any advice would be welcome, thank you.

2

Since A is real, x^TAx = x^T \frac{(A+A^T)}{2}x. Thus, since x^TAx/x^TBx is a generalized Rayleigh quotient, this should be equivalent to a generalized eigenvalue problem involving the symmetric part of A (up to the normalization constraint).

4 Likes
3

In particular, the maximum occurs when x is an eigenvector corresponding to the largest eigenvalue λ of the real-symmetric generalized eigenproblem (A+A^T)x=2Bx.

You can solve this by using eigen(Symmetric(A+A'), Symmetric(2B)) in Julia if the matrices are not too large and looking at the last eigenvalue and eigenvector (since the eigenvalues are sorted), or using (e.g.) eigs in ARPACK.jl for large sparse matrices.

5 Likes
4

This is a generalized Rayleigh quotient. It is normally used for generalized eigenvalue problems (in this case Ax = \lambda Bx) but that requires A to be symmetric. The expression x^\top A x only depends on the symmetric part of A. Nonsymmetric eigenvalue problems are much more difficult than symmetric ones, you could look at some references on that. For the case where A^\top = A, an Arnoldi method is a popular choice. This is not an optimization method but some type of power method for eigenvalue problems.

1 Like
5

Thank you for your replies! I was not expecting a solution as elegant as this.

I coded a minimal example:

using LinearAlgebra

n = 4
A = randn(n, n)
rootB = randn(n, n)
B = rootB'rootB

eig = eigen(Symmetric(A+A'), Symmetric(2B))
x = eig.vectors[:, n]
x'x

I notice the generalized eigenvectors don’t have unit norm. Can I simply scale them so this constraint is satisfied?

6

Welcome to the Discourse community :slight_smile:.

Yes - you can scale eigenvectors x = x / (x'*x) to satisfy x^Tx=1, as this shouldn’t change the generalized Rayleigh quotient. For a general scaling, y = \alpha x, and \frac{y^TAy}{y^TBy} = \frac{\alpha^2 x^TAx }{\alpha^2 x^TBx} = \frac{x^TAx }{x^TBx}.

1 Like