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

share|improve this question
1  
Use _?NumericQ in your function definition –  rm -rf May 13 '13 at 20:10

1 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 ---

share|improve this answer
    
Thank you for such a great answer! One thing though: the pnorm you have here outputs a different answer than my original one. It seems like yours just raises the elements to the pth power. Doesn't a p-norm also involve taking the pth root of the sum? –  Muhammad Khan May 13 '13 at 21:34
    
Yes it does. Minimizing one is equivalent to minimizing the other (it's a monotonic function) but this change of mine messes up the result. Will edit accordingly. –  Daniel Lichtblau May 13 '13 at 22:23
    
It seems quite unstable and can give hugely negative results when I change to account for pth roots. Is there any chance that the additive term was meant to be lambda*Tr[mm.c.(Transpose[mm])]? That seems to make it give plausible results quickly. –  Daniel Lichtblau May 13 '13 at 22:56
    
Well the model I'm trying to minimize is sum_i(sum_j(c_ij * mm_ij)) <- this is equivalent to Tr[c.(Transpose[mm])] But I do see what you mean about the pth root being unnecessary if I only want argmin –  Muhammad Khan May 13 '13 at 23:21
    
Okay. I'm still getting results that do not seem to be playing nicely with the code. Will have another look tomorrow. –  Daniel Lichtblau May 13 '13 at 23:32

Your Answer

 
discard

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.