Beefy Boxes and Bandwidth Generously Provided by pair Networks
Don't ask to ask, just ask

Help with Matrix math!

by madslynge
on Dec 16, 2007 at 18:10 UTC ( #657293=perlquestion: print w/replies, xml ) Need Help??
madslynge has asked for the wisdom of the Perl Monks concerning the following question:

Hey i need some help with maximizing a matrix
max (x) x'Dx
st. x'x=1
D is a symmetrical (n x n)matrix and x is an unknown (n x 1) vector

Replies are listed 'Best First'.
Re: Help with Matrix math!
by ikegami (Pope) on Dec 16, 2007 at 18:56 UTC
    You could use PDL, an efficient math library for Perl. It can do matrix math with ease.
Re: Help with Matrix math!
by ww (Archbishop) on Dec 16, 2007 at 18:24 UTC
    What, if anything, have you tried?

    Is there a Perl component to this question or did you mistake the Monastery for the Cloistered Towers of MathMonks?

    and, please, don't double-post your questions.

Re: Help with Matrix math!
by pc88mxer (Vicar) on Dec 16, 2007 at 21:16 UTC
    This sounds a lot like a homework problem, but I'll give you the following pointers:
    1. This is equivalent to finding the largest eigenvalue of D.
    2. A standard way of solving this problem is to use the "power method" which is documented here: wikipedia entry
    3. The perl module PDL is useful for performing matrix operations.
Re: Help with Matrix math!
by randyk (Parson) on Dec 16, 2007 at 21:33 UTC

    If you set aside for a moment the constraint on the vector x, then the problem of maximizing xiDijxj (a sum over repeated indices is implied) translates, upon taking the derivative with respect to xi, into the equation Dijxj = 0. This is a special form of the eigenvalue equation for the matrix D, with the eigenvalue being 0. In order for non-trivial solutions, the matrix D cannot have an inverse, so it's determinant must vanish. Finding the eigenvectors can be done with PDL; see Mastering Algorithms with Perl for a discussion. The Math::Cephes::Matrix module can also do this for real symmetric matrices. Generally, the components of the eigenvectors found in this way are not completely determined; by convention, the normalization xixi = 1 is imposed.

    Update: The preceding needs to be augmented by a test on the Hessian of the matrix to determine what type of critical point is present.

      In this case you have to take into account the constraint on the vector x. The quadratic form xiDijxj is unbounded on Rn, so a derivative test (even with an Hessian test) is only going to find local extreme points, and there is no guarantee that those points will lie on the surface of interest (i.e. x'x = 1.) The difference is that for a point to be a local maximum when x is constrained, the derivative only needs to be zero when evaluated in the tangent space of the constraining surface at that point, not zero in every direction.

      The standard way to include the surface constraint is to use Lagrangian multipliers.

        I was thinking of the problem differently - find the extremal points of xDx, which leads to the eigenvalue equation Dx = 0. This equation doesn't determine xx completely, so one is free to impose, for example, xx = 1 as a normalization condition. But you're right that if xx = 1 is intended as a true constraint, then a method like Lagrange multipliers should be used.
Re: Help with Matrix math!
by CountZero (Bishop) on Dec 16, 2007 at 20:12 UTC
    How would you do it by hand / paper & pencil? Do you apply a certain algorithm? If so, that algorithm can then most probably be translated into a computer program.


    A program should be light and agile, its subroutines connected like a string of pearls. The spirit and intent of the program should be retained throughout. There should be neither too little or too much, neither needless loops nor useless variables, neither lack of structure nor overwhelming rigidity." - The Tao of Programming, 4.1 - Geoffrey James

Re: Help with Matrix math!
by rg0now (Chaplain) on Dec 17, 2007 at 19:04 UTC
    What you are trying to solve is a nonlinearly constrained nonlinear (actually, quadratic) mathematical program. This is a tough one, I can tell you...

    However, all is not lost, because you might be able to simplify your problem into something actually soluble.

    First: are you absolutely sure that you want to optimize over the boundary of the unit ball (this is what the constraint x'x=1 seems to impose)? Isn't it enough to only assure that x is constrained into a sane region, like e.g., the n-dimensional unit box? If it is, you get a quadratic program:

    max_x (x' D x) s.t. -1 <= x_i <= 1 \forall i=1..n
    See more on quadratic programs here: Quadratic programming.

    Second: all depends on whether or not your D is negative definite or not (all eigenvalues are negative). If not, you are out of luck: the problem is (at least) NP hard. If, on the other hand, D is nicely negative definite, then the objective function is concave and the uniqueness of the optimal solution is guaranteed. In this case, you could use a commercial quadratic solver, like CPLEX, or, after suitable linearization, a linear program solver. too. E.g., this one Math::GLPK has a nice Perl interface (believe me it's nice: I wrote it!...:-).

Log In?

What's my password?
Create A New User
Node Status?
node history
Node Type: perlquestion [id://657293]
Approved by ww
and all is quiet...

How do I use this? | Other CB clients
Other Users?
Others cooling their heels in the Monastery: (5)
As of 2018-05-25 02:03 GMT
Find Nodes?
    Voting Booth?