# Using FindMinimum with matrix input

I define the following function which spits out a number S when I feed it a dxd square matrix:

energycoeffs[coeffmatrix_] :=
Module[{psi, psinorm, U, P, dxU, dyU, L, S},
psi[x_, y_] := coeffmatrix.thetabasis[x, y];
psinorm[x_, y_] := Chop[Norm[psi[x, y]]];
U[x_, y_] := 1/psinorm[x, y]* psi[x, y];
P[x_, y_] := U[x, y].ConjugateTranspose[U[x, y]];
dxU[x_, y_] := (U[x + dx/2, y] - U[x - dx/2, y])/dx /. dx -> 0.01 ;
dyU[x_, y_] := (U[x, y + dy/2] - U[x, y - dy/2])/dy /. dy -> 0.01 ;
L[x_, y_] :=
Re[ConjugateTranspose[dxU[x, y]].dxU[x, y] +
ConjugateTranspose[dyU[x, y]].dyU[x, y] -
ConjugateTranspose[dxU[x, y]].P[x, y].dxU[x, y] -
ConjugateTranspose[dyU[x, y]].P[x, y].dyU[x, y]];
S :=
N1*N2*NIntegrate[
L[x, y]*Boole[{x, y} ∈ ℛ], {x, 0, gam1 + Re[gam2]}, {y, 0, Im[gam2]},
Method -> "AdaptiveMonteCarlo", AccuracyGoal -> 3];
S]


I evaluate it a point and it works:

energycoeffs[{{1, 0, 0}, {0, 0, 0}, {0, 0, 1}}]

{{18.8666}}


But I try to use it to find the coefficients of the optimal diagonal matrix,

FindMinimum[
{
energycoeffs[{{p, 0, 0}, {0, q, 0}, {0, 0, r}}],
0 <= p <= 1.2, 0 <= q <= 1.2, 0 <= r <= 1.2
},
{p, q, r}]


it gives me the following error over and over again:

NIntegrate::inumr: The integrand Boole[{x, y} ∈ Polygon[{{0, 0}, {0, Power[<<2>>] π}, {Power[<<2>>] π, Power[<<2>>] π}, {Power[<<2>>] π, 0}}]] Re[<<1>>] has evaluated to non-numerical values for all sampling points in the region with boundaries {{0, 1}, {0, 1}}. >>

Does anyone know what I am doing wrong here?

EDIT: Here is the full code leading up to the function.

(* Set parameters *)
d = 3; (* !! *)
N1 = 1; N2 = 1; (*Number of unit \
cells*)
tau =
Exp[I π/3]; (* γ2=tau*γ1*)
gam1 = π Sqrt[
d]; gam2 = π Sqrt[d] tau;
thetap[x_, y_, p_] = (
E^(-((I (Sqrt[d] π tau - 2 (x + I*y))^2)/(4 π tau)))
EllipticTheta[
3, -(((d + 2 p) π)/(2 d)) + (x + I*y)/(Sqrt[d] tau),
E^(-((I π)/(d tau)))])/Sqrt[-I d tau];
thetabasis[x_, y_] = Table[{thetap[x, y, p]}, {p, 0, d - 1}];
Dimensions[thetabasis[x, y]]

ℛ =
Polygon[{{0, 0}, {Re[gam2], Im[gam2]}, {gam1 + Re[gam2],
Im[gam2]}, {gam1, 0}}] ;
Graphics[ℛ]


UPDATE: c186282 hit the nail on the head. The output of the function was not a scalar but a list with the scalar as the [1,1] value. E.g.-

energycoeffs[{{1, 0, 0}, {0, 0, 0}, {0, 0, 1}}]

{{18.8666}}


So now, I just use Norm[] to extract the scalar value. However, now I have a further question. FindMinimum is no longer throwing up any error messages, but it is extremely slow! It's been running for 10 minutes now without giving an answer. Would anyone have any advice as to how to fix it? Is it a problem with how I define the function?

UPDATE: I tried out the solution suggested by @george2079, and it works pretty well. I created a discretised mesh.

nsamp = 14.;
dx = (gam1 + Re[gam2])/nsamp  ; dy =
Im[gam2]/nsamp ;(* Equivalent dx value. Note: x and y directions are \
not same, depending on shape of the unit cell. *)

axis = Flatten[Table[{i*dx, j*dy}, {i, 0, nsamp}, {j, 0, nsamp}],
1] ;
Dimensions[axis]

(* Select points in unit cell, always <= L1 *)
pts = Select[axis, # \[Element] \[ScriptCapitalR] &];
Dimensions[pts]
ListPlot[pts]


Then I copy and pasted his newly defined Module[], with a slight tweak in the integration definition:

energycoeffs[coeffmatrix_] :=
Module[{psi, psinorm, U, P, dxU, dyU, L, S},
psi[x_, y_] := coeffmatrix.thetabasis[x, y];
psinorm[x_, y_] := Chop[Norm[psi[x, y]]];
U[x_, y_] := 1/psinorm[x, y]*psi[x, y];
P[x_, y_] := U[x, y].ConjugateTranspose[U[x, y]];
dxU[x_, y_] :=
With[{dx = 0.01}, (U[x + dx/2, y] - U[x - dx/2, y])/dx];
dyU[x_, y_] :=
With[{dy = 0.01}, (U[x, y + dy/2] - U[x, y - dy/2])/dy];
L[x_, y_] :=
First@First@
Re[ConjugateTranspose[dxU[x, y]].dxU[x, y] +
ConjugateTranspose[dyU[x, y]].dyU[x, y] -
ConjugateTranspose[dxU[x, y]].P[x, y].dxU[x, y] -
ConjugateTranspose[dyU[x, y]].P[x, y].dyU[x, y]];
N1 N2 Total@(L @@@ pts) * dx *dy]


It works pretty well!

In[80]:= energycoeffs[IdentityMatrix[d]]  // Timing

Out[80]= {7.073925, 19.6902}


Interestingly, it takes 7x more time to get approximately the same answer that we got with randomly generated points!

• can you include the definition for thetabasis, gam1, gam2, N1, and N2 so I could give it a try? – Jason B. Feb 9 '15 at 13:51
• I think you should be looking at what happens when energycoeffs[{{p, 0, 0}, {0, q, 0}, {0, 0, r}}] is evaluated. – m_goldberg Feb 9 '15 at 14:04
• If I scrape and paste your example directly then one problem I am seeing is Boole[{1, 2} \[Element] \[ScriptCapitalR]]. Mathematica does not understand that you dearly want Olde English reverse italic animated font with lots of cute curlies "R" to mean "Reals". When I replace that with Boole[{1, 2} \[Element] Reals] Mathematica understands you are asking if 1 and 2 are elements of the real numbers. I hope there is someone inside Wolfram who cries on a daily basis as a result of new typesetting features that make users think they MUST desktop publish EVERY detail of their Mathematica input. – Bill Feb 9 '15 at 19:47
• Obviously I cannot run your code either but here are some ideas. Make sure that your function is returning a scalar and force your inputs to be numbers. energycoeffs[p_?NumericQ,q_?NumericQ,r_?NumericQ]:=energycoeffs[{{p, 0, 0}, {0, q, 0}, {0, 0, r}}] Then minimize energycoeffs[p,q,r] – c186282 Feb 10 '15 at 1:02
• I have edited the main post, but I have a side question as well. How exactly does the 'p_?NumericQ' syntax work? Is it some kind of typecasting? – ap21 Feb 11 '15 at 19:47

Here is an approach that is much faster, but with potentially some loss of accuracy. Do away with NIntegrate and eval the function at preselected points in the region:

 nsamp = 200;
pts = Select[
Table[{RandomReal[{0, gam1 + Re[gam2]}],
RandomReal[{0, Im[gam2]}]}, {nsamp}] , # \[Element] \[ScriptCapitalR] & ];
energycoeffs[coeffmatrix_] :=
Module[{psi, psinorm, U, P, dxU, dyU, L, S},
psi[x_, y_] := coeffmatrix.thetabasis[x, y];
psinorm[x_, y_] := Chop[Norm[psi[x, y]]];
U[x_, y_] := 1/psinorm[x, y]*psi[x, y];
P[x_, y_] := U[x, y].ConjugateTranspose[U[x, y]];
dxU[x_, y_] := With[{dx = 0.01}, (U[x + dx/2, y] - U[x - dx/2, y])/dx];
dyU[x_, y_] := With[{dy = 0.01}, (U[x, y + dy/2] - U[x, y - dy/2])/dy];
L[x_, y_] := First@First@
Re[ConjugateTranspose[dxU[x, y]].dxU[x, y] +
ConjugateTranspose[dyU[x, y]].dyU[x, y] -
ConjugateTranspose[dxU[x, y]].P[x, y].dxU[x, y] -
ConjugateTranspose[dyU[x, y]].P[x, y].dyU[x, y]];
N1 N2 Total@(L @@@ pts) (gam1 + Re[gam2]) Im[gam2] / nsamp ]
energycoeffs[{{1, 0, 0}, {0, 0, 0}, {0, 0, 1}}] // Timing


{1.263608, 19.7669}

compare with ,

{24.429757, 18.6672}

a 20 fold performance increase. I don't know how significant that error is for your application, but you might use this to get close to the minimum then use the integral form for refinement. ( The accuracy of this approach will be improved greatly if you can discretize (mesh) your region, instead of using random points )