Minimizing a Matrix

Question

pnorm[A_, p_] := Module[{m, n, x, y, f, ans, ret},
  {m, n} = Dimensions[A];
  y = Array[x, n];
  f = Norm[A.y, p]/Norm[y, p] // FullSimplify;
  ans = Maximize[f, Array[x, n]] // N;
  ret = Part[ans, 1];
  ret
]

findD[L_, c_, \[Lambda]_, p_] := Module[{B, M, w, pp, ans},
  M = Array[B, Dimensions[L]];
  w = Tr[c.(Transpose[M])];
  pp = pnorm[L - M, p] // FullSimplify;
  ans = Minimize[Hold[pp + \[Lambda] w], 
    Flatten[Array[B, Dimensions[L]]]] // N;
  ans
]

findD[{{3, 8}, {2, 5}}, {{9, 1}, {3, 6}}, .5, 3]

During evaluation of In[44]:= NMinimize::nnum: The function value 6.59387       +1/((Abs[<<1>>]^3+<<1>>^3)/(Abs[<<1>>]^3+<<1>>^3))^(1/3) is not a number at {B$33209[1,1],B$33209[1,2],B$33209[2,1],B$33209[2,2]} = {0.643476,0.766972,0.280303,0.964763}. >>

During evaluation of In[44]:= NMinimize::nnum: The function value 1/((Abs[<<1>>]^3+<<1>>^3)/(Abs[<<1>>]^3+<<1>>^3))^(1/3)+0.5 (9 B$33209[1,1]+B$33209[1,2]+3 B$33209[2,1]+6 B$33209[2,2]) is not a number at {B$33209[1.,1.],B$33209[1.,2.],B$33209[2.,1.],B$33209[2.,2.]} = {0.643476,0.766972,0.280303,0.964763}. >>

During evaluation of In[44]:= NMinimize::nnum: The function value 1/((Abs[<<1>>]^3+<<1>>^3)/(Abs[<<1>>]^3+<<1>>^3))^(1/3)+0.5 (9 B$33209[1,1]+B$33209[1,2]+3 B$33209[2,1]+6 B$33209[2,2]) is not a number at {B$33209[1.,1.],B$33209[1.,2.],B$33209[2.,1.],B$33209[2.,2.]} = {0.643476,0.766972,0.280303,0.964763}. >>

During evaluation of In[44]:= General::stop: Further output of NMinimize::nnum will be suppressed during this calculation. >>

Minimize[
 Hold[pp$33209 + 0.5 w$33209], {B$33209[1., 1.], B$33209[1., 2.], 
 B$33209[2., 1.], B$33209[2., 2.]}]

The code above is trying to find a matrix $\mathbf M$ which minimizes $\|\mathbf L-\mathbf M\|_p + \lambda\,\mathrm{tr}(\mathbf C\mathbf M^\top)$. However, the usage of Minimize in the findD function is giving me problems. Does anyone know how to fix the "not a number" problem?

Thank you!

Answer 1

You will need to make pnorm into a black-box function so that it never tries to evaluate symbolically. Also it is probably worthwhile to use NMinimize/NMaximize in both. Finally the minimax process seems to behave better if you take pth powers in the inner optimization.

pnorm[aa_, p_] /; MatrixQ[aa, Element[N[#], Reals] &] := Module[
  {m, n, x, y, f, ans},
  {m, n} = Dimensions[aa];
  y = Array[x, n];
  f = Total[Abs[aa.y]^p]/Total[Abs[y]^p];
  ans = NMaximize[f, y];
  ans[[1]]]

findD[ll_, c_, \[Lambda]_, p_] := Module[
  {bb, mm, w, ans},
  mm = Array[bb, Dimensions[ll]];
  w = Tr[c.(Transpose[mm])];
  pp = pnorm[ll - mm, p];
  ans = NMinimize[pnorm[ll - mm, p] + \[Lambda] w, Flatten[mm]];
  ans]

Your example:

findD[{{3, 8}, {2, 5}}, {{9, 1}, {3, 6}}, .5, 3]

(* {86.9382369925, {bb$7195688[1, 1] -> 2.99883475189, 
      bb$7195688[1, 2] -> 4.10214829344, 
  bb$7195688[2, 1] -> -1.89824492827, 
      bb$7195688[2, 2] -> 5.00101503916}} *)

--- edit ---

Here is code that seems to work better.

pnorm[aa_, p_] /; MatrixQ[aa, Element[N[#], Reals] &] := Module[
  {m, n, x, y, f, ans},
  {m, n} = Dimensions[aa];
  y = Array[x, n];
  f = Total[Abs[aa.y]^p]/Total[Abs[y]^p];
  ans = FindMaximum[f, y];
  {ans[[1]]^(1/p), y /. ans[[2]]}]

findD[ll_, c_, \[Lambda]_, p_] := Module[
  {bb, mm, obj, ans},
  mm = Array[bb, Dimensions[ll]];
  obj[mat_?(MatrixQ[#, Element[N[#], Reals] &] &), lmat_, cmat_, 
    pval_, lam_] :=
   Module[{res = pnorm[lmat - mat, pval]},
    res[[1]] + lam*Total[Abs[Flatten[cmat*Transpose[mat]]]]];
  ans = FindMinimum[Evaluate[obj[mm, ll, c, p, \[Lambda]]], 
    Flatten[mm]];
  ans]

findD[{{3, 8}, {2, 5}}, {{9, 1}, {3, 6}}, .5, 3]

(* Out[201]= {9.93991501159, {bb$39780695[1, 1] -> 0.00147967770299, 
      bb$39780695[1, 2] -> 0.0190392337428, 
  bb$39780695[2, 1] -> -2.33558050383*10^-9, 
      bb$39780695[2, 2] -> -0.00620592593044}} *)

One thing to note is that I used absolute values for terms in the second part of the objective function. Else i think it can be arbitrarily small.

I also return the maximizing values for the inner optimization. This is from prior experiments, not shown, that also figured sizes of vectors in the objective, to penalize for getting too small. That seems not to be needed in this most recent version but thought I'd leave that there in case you do decide to account for such vector sizes.

--- end edit ---