Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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_] := 
    (* 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 *)
      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?

share|improve this question
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. –  Rahul Narain May 22 '13 at 1:21
show 4 more comments

2 Answers

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

share|improve this answer
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
add comment

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.

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

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

share|improve this answer
add comment

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.