Optimizing functions taking matrix arguments

Question

I'm looking for general information on how to optimize matrix valued functions, I have the following function I'm looking to maximize (or figure out if this is possible at all).

MaximizeFunction[W_, DataCoupled_] := 
  Module[{newDataCouple},
    (* Elementwise Multiplication on a list of 2 element vectors *)
    newDataCouple = Flatten[List[Dot[W, #] & /@ DataCoupled], 1];
    (* Take the first and second elements of each of the vectors in the previous list and 
      perform an independence test on them to obtain the p-value *)
    Return[
      IndependenceTest[Extract[#, 1] & /@ newDataCouple,Extract[#, 2] & /@ newDataCouple]]];

Input values are of the type:

DataCouple = {{.5, .8}, {.7, .9}, {.6, .9}, ... } 
W = {{1, 0}, {0, 1}}

Could I then use NMaximize or Maximize to optimize this function?

What are the free variables? Are you trying to find the W such that this function is maximized? Are there any constraints on W? — Jonathan Shock, May 22 '13 at 00:57
If it's a matter of finding W such that this function is closest to 1, there are a whole family of solutions. You can define a function tomax[w1_?NumberQ, w2_?NumberQ, w3_?NumberQ, w4_?NumberQ] := MaximizeFunction[{{w1, w2}, {w3, w4}}, DataCouple] and use NMaximize on this. However, plotting with any values of w2,w3,w4 you will find that there will always be a value of w1 which maximizes the function (to 1). (This is very inelegant coding but it should work). — Jonathan Shock, May 22 '13 at 01:04
No constraints on W, but I know that the output is invariant under scaling of W. One proposed solution to this is to constrain W such that WW*= 2x2 Identity. The number of elements in DataCouple are high (~120000), but I will try what you suggested. — Ryan Warnick, May 22 '13 at 01:08
Would you be happy with a single solution, or do you want the whole family? If you just want a single solution then it appears from this example to be an optimization problem in one variable. — Jonathan Shock, May 22 '13 at 01:10
This isn't a matrix-valued function. It's a function whose arguments are matrices, but the value the function returns is a scalar $p$-value. — , May 22 '13 at 01:21
Return is not needed here and you can also drop the ;. This is not just nitpicking: while in this particular example it does not make a difference if you use it or not, when used inside some functions, it will return from those functions and not from your function. — Szabolcs, May 22 '13 at 01:55
This is not an answer to your question, but I believe you can replace: Flatten[List[Dot[W, #] & /@ DataCoupled], 1] with: DataCoupled.W — Mr.Wizard, May 22 '13 at 06:25

Jonathan Shock · Answer 1 · 2013-05-22T01:40:16.143

3

(Edit, I've edited the following almost entirely from the original, but the idea remains the same)

From the comments it seems that a single solution will be enough. You want the input of the original function to be a numerical matrix. You can set up a test for this as follows:

matrixnumQ[exp_] := MatrixQ[exp, NumericQ]

Then defining your original function to test this on the input you can use the original function directly in NMaximize.

MaximizeFunction[W_?matrixnumQ, DataCoupled_] := Module[{newDataCouple},
newDataCouple = Flatten[List[Dot[W, #] & /@ DataCoupled], 1];
Return[IndependenceTest[Extract[#, 1] & /@ newDataCouple, 
 Extract[#, 2] & /@ newDataCouple]]];

DataCouple = {{.5, .8}, {.7, .9}, {.6, .9}, {.4, 0.8}}
NMaximize[MaximizeFunction[{{w1, 1}, {1, 1}}, DataCouple], {w1}]

This should give you an answer without any errors.

edited May 22 '13 at 01:40

answered May 22 '13 at 01:19

Jonathan Shock

3,015
15
24

When I put in W_?NumberQ for the MaximizeFunction I receive the following error immediately: The function value MaximizeFunction[...] is not a number at {w1}=.0914636 – Ryan Warnick May 22 '13 at 01:35
@RyanWarnick I've edited the answer to take this into account. This should solve the problem, but I'm sure that there is a more elegant way to implement the test of whether the input is a number-valued matrix. – Jonathan Shock May 22 '13 at 01:41

Michael E2 · Answer 2 · 2013-05-23T01:07:28.813

See if this is what you're after. I took the liberty of simplifying your MaximizeFunction, and in the process it became about twice as fast. I also got rid of the initial capitals. Best to avoid them, and avoid conflicting inadvertently with built-in functions.

In a comment you indicate that it might be sufficient to find the maximum over orthogonal matrices ($WW^* = I$). Then it is easy to do, since RotationMatrix[t] gives half of them. The other half of the orthogonal matrices are given by any reflection times a rotation matrix, such as {{1, 0}, {0, -1}}.RotationMatrix[t]. In all cases I tried, the p-value, as a function of t was the same for reflections as for rotations; further, the period as a function of t was π/2. (Perhaps one should check that.) If so, we can just use rotations.

maximizeFunction[W_, DataCoupled_] := Module[{newDataCouple},
   newDataCouple = DataCoupled.Transpose[W];
   IndependenceTest[First /@ newDataCouple, Last /@ newDataCouple]];
obj[t_?NumericQ, couple_] := maximizeFunction[RotationMatrix[t], couple]

We'll make up a large, random data set. It turns out there can be several local maxima, so using FindMaximum would probably give unreliable results. Another problem is that it takes a long time to evaluate a single function call. This makes using NMaximize take a very long time.

SeedRandom[1];
dc2 = RandomReal[{0, 1}, {120000, 2}]
(Table[obj[t, dc2], {t, 0.1, 1., 0.1}]; // AbsoluteTiming // First) / 10

0.2092696

Plot the function to get a sense of where the maximum is. (Plot from 0 to 2 Pi to check the periodicity.)

plot = Plot[obj[t, dc2], {t, 0, \[Pi]/2}, MaxRecursion -> 1]

One period of obj[t, dc2]

We can get a rough approximation of the maximum from plot:

maxpt = Last@SortBy[Cases[plot, {_Real, _Real}, Infinity], Last]

{0.897961, 0.993164}

Use the first coordinate as an initial point for FindMaximum.

t0 = maxpt[[1]];
({pvalue, tsol} = FindMaximum[obj[t, dc2], {t, t0, t0 + 1/100}]) // AbsoluteTiming

{7.493610, {1., {t -> 0.891873}}}

Since we got a p-value of 1, we know it's the maximum. Here is the optimal $W$:

RotationMatrix[t] /. tsol

{{0.627955, -0.778249}, {0.778249, 0.627955}}

Here's a function that does the whole thing:

findMax[couple_] := Block[{plot, t0},
  plot = Plot[obj[t, couple], {t, 0, \[Pi]/2}, MaxRecursion -> 1];
  t0 = First @ Last @ SortBy[Cases[plot, {_Real, _Real}, Infinity], Last];
  FindMaximum[obj[t, couple], {t, t0, t0 + 1/100}]]

findMax[dc2] // AbsoluteTiming

{26.854463, {1., {t -> 0.891873}}}

If maximizing over rotation matrices is not sufficient, then you could do something similar with a different parametrization of the matrices. It tends to get harder as the dimension of the input domain increases.

If you know a formula for the p-value of the IndependenceTest, you might be able to use that to speed things up. (If there is a formula that can be differentiated, then FindMaximum can use Newton's method and so on.)

Optimizing functions taking matrix arguments

2 Answers2

Linked