# Optimizing functions taking matrix arguments

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 0: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 1: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 1: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 1: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. – user484 May 22 '13 at 1:21

(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.

• 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 1: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 1:41

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;
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] 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[];
({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.)