I am trying to solve a (Sequence) Quadratic Programming (SQP) problem in Mathematica.

Here’s the problem setup:

size = 9; (* for illustrative purposes, can be larger *)
(* known coefficients *)
m1 = {{a, b, c, d}, {e, f, g, h}, {i, j, k, l}};
m2 = {{m, 0, 0, 0}, {0, n, 0, 0}, {0, 0, o, 0}};
(* unknowns *)
θ = Table[{Subscript[x, i], Subscript[y, i], Subscript[z, i], 1}, {i, size}];

(* My system of quadratic equations *)
data == Transpose [m1.Transpose[θ] + m2.Transpose[θ^2]];

(*Subject to the following constraints *)
constraints =
Table[{Subscript[x, i]^2 + Subscript[y, i]^2 + Subscript[z, i]^2 == 1}, {i, size}];


Sample data:

data =
{{0.0326857, -0.23977, -0.980886}, {-0.00491715, 0.160488, -0.980237},
{-0.0189063, 0.137756, 0.969873}, {-0.0165821, -0.256267, 0.967219},
{-0.00898421, -1.00759, -0.258025}, {-0.0131764, 0.918642, 0.192466},
{-0.00673595, 0.0658302, -0.996648}, {0.00647973, -0.147268, -0.997363},
{-0.0200075, -0.164866, 0.981011}, {-0.0242359,0.0523168, 0.981269},
{-0.0132569, -1.03757, -0.0196142}, {-0.0164502, 0.94126, -0.0110496},
{-1.01661, -0.0493139, -0.0169287}, {-0.0250339, -0.0569752, 0.988069},
{0.981252, -0.048789, -0.00833748}, {-0.00965323,-0.0403172, -1.00421}};

m1 =
{{0.851943, 0, 0, 0.0154734}, {0.500424, 1.01131, 0, 0.0564313},
{-0.164399, 0.00390119, 1.00309, 0.00583693}};

m2 = {{0.00912, 0, 0, 0}, {0, -0.00320, 0, 0}, {0, 0, 0.00175, 0}};

• Added some sample data to help...
– Pam
Commented Sep 9, 2013 at 22:06
• What is your objective function? What you have above is a matrix, not a function. As written, it will give a recursion error, but that's a different matter: should not reuse f in defining it since f is an element of m2. Commented Sep 9, 2013 at 22:41
• Daniel, I have a set of quadratic equations given by data == Transpose [m1.Transpose[θ] + m2.Transpose[θ^2]];
– Pam
Commented Sep 9, 2013 at 23:23
• Edited the original question (f is not re-used any more).
– Pam
Commented Sep 9, 2013 at 23:26
• @Daniel: There is not a single objective function but a set of equations subject to a set of constraints...
– Pam
Commented Sep 9, 2013 at 23:28

In your problem we can try matrix norm minimization. One straight forward way will be to use FindMinimum as follows

size = 16;
θ =Transpose[(Array[#, size] & /@ {x, y, z})~Join~{ConstantArray[1, size]}];
cons = Flatten[θ /. {a_, b_, c_,1} -> {a^2 + b^2 + c^2 == 1}] /. List :> And;
res = FindMinimum[{
Norm[Transpose[m1.Transpose[θ] +m2.Transpose[[Theta]^2]]-data,"Frobenius"],
cons},
Transpose[{Flatten@(θ[[All, 1 ;; 3]]),RandomReal[{-1, 1}, 3 size]}]
];


The local minimum is not too bad

First@res


0.374061

do a recheck as

ob=Transpose[m1.Transpose[θ]+m2.Transpose[θ^2]]-data;
Norm[ob/.res[[2]],"Frobenius"]


0.374061

Now we can see that constraints are fulfilled with a good tolerance

conResolve = (θ /. res[[2]] /. {a_, b_, c_, 1} -> a^2 + b^2 + c^2 - 1);
{Mean@conResolve, Variance@conResolve}


{3.38641*10^-13,1.26468*10^-24}

With a random assignment of solution your given objective tend to behave like the following

LaunchKernels[10];
SeedRandom[1234];
check = BlockRandom[
vals = ParallelTable[RandomReal[{-10, 10}, 2], {10^5}];
ParallelMap[(Flatten[{#,Norm[ob /. (Transpose[{Flatten@(θ[[All, 1 ;; 3]]),
RandomReal[Sort@#, 3 size]}] /. {a_, b_} :> Rule[a, b]),"Frobenius"]
}] &) /@ # &, Partition[vals, 10^4]]];


The minimum value of the objective I could get with these random solution assignment from different random sub-intervals of $[-10,10]$.

Last@First@SortBy[Flatten[check, 1], Last]


3.80535

is pretty bad compared to what FindMinimum gave us out of the box.

You can see the acceptable trend for the occurrence of a local minimum which looks like a global one for the considered interval! This is bad and can stall any optimization algorithm.

• Brilliant. This is more than I needed.
– Pam
Commented Sep 10, 2013 at 13:11
• Kudos to Daniel as well for figuring out that f was an element of m2 and I was repeating a variable in two instances.
– Pam
Commented Sep 10, 2013 at 13:12

Correcting a few things, it can be set up as follows.

size = 16;
\[Theta] =
Table[{Subscript[x, i], Subscript[y, i], Subscript[z, i], 1}, {i,
size}];
constraints =
Table[{Subscript[x, i]^2 + Subscript[y, i]^2 + Subscript[z, i]^2 ==
1}, {i, size}];

data = {{0.0326857, -0.23977, -0.980886}, {-0.00491715,
0.160488, -0.980237}, {-0.0189063, 0.137756,
0.969873}, {-0.0165821, -0.256267,
0.967219}, {-0.00898421, -1.00759, -0.258025}, {-0.0131764,
0.918642, 0.192466}, {-0.00673595,
0.0658302, -0.996648}, {0.00647973, -0.147268, -0.997363}, \
{-0.0200075, -0.164866, 0.981011}, {-0.0242359, 0.0523168,
0.981269}, {-0.0132569, -1.03757, -0.0196142}, {-0.0164502,
0.94126, -0.0110496}, {-1.01661, -0.0493139, -0.0169287}, \
{-0.0250339, -0.0569752,
0.988069}, {0.981252, -0.048789, -0.00833748}, {-0.00965323, \
-0.0403172, -1.00421}};

m1 = {{0.851943, 0, 0, 0.0154734}, {0.500424, 1.01131, 0,
0.0564313}, {-0.164399, 0.00390119, 1.00309, 0.00583693}};

m2 = {{0.00912, 0, 0, 0}, {0, -0.00320, 0, 0}, {0, 0, 0.00175, 0}};


Now define the polynomials we want to approximately satisfy.

polys = data -
Transpose[m1.Transpose[\[Theta]] + m2.Transpose[\[Theta]^2]];

fpolys = Flatten[polys];


To optimize, we could minimize the sum of squares of fpolys. I do that below. Unfortunately this is no longer a quadratic objective function, but it is still amenable to e.g. nonlinear interior point methods (the presence of constraints rule out some other approaches).

{min, vals} = FindMinimum[{fpolys.fpolys, Flatten[constraints]},
Evaluate[Sequence @@ Variables[polys]]]

(* Out[16]= {0.139921, {Subscript[x, 1] -> 0.0199536,
Subscript[x, 2] -> -0.0251878, Subscript[x, 3] -> -0.0358971,
Subscript[x, 4] -> -0.0362198, Subscript[x, 5] -> -0.0588463,
Subscript[x, 6] -> -0.0943821, Subscript[x, 7] -> -0.0244723,
Subscript[x, 8] -> -0.00827959, Subscript[x, 9] -> -0.0374074,
Subscript[x, 10] -> -0.0409569, Subscript[x, 11] -> -0.063582,
Subscript[x, 12] -> -0.0958467, Subscript[x, 13] -> -0.939709,
Subscript[x, 14] -> -0.0419559, Subscript[x, 15] -> 0.878312,
Subscript[x, 16] -> -0.0260752, Subscript[y, 1] -> -0.294872,
Subscript[y, 2] -> 0.115827, Subscript[y, 3] -> 0.102224,
Subscript[y, 4] -> -0.292127, Subscript[y, 5] -> -0.965338,
Subscript[y, 6] -> 0.978937, Subscript[y, 7] -> 0.0212626,
Subscript[y, 8] -> -0.193187, Subscript[y, 9] -> -0.202956,
Subscript[y, 10] -> 0.0166422, Subscript[y, 11] -> -0.997514,
Subscript[y, 12] -> 0.994625, Subscript[y, 13] -> 0.306808,
Subscript[y, 14] -> -0.0936467, Subscript[y, 15] -> -0.464459,
Subscript[y, 16] -> -0.0814992, Subscript[z, 1] -> -0.955328,
Subscript[z, 2] -> -0.99295, Subscript[z, 3] -> 0.994114,
Subscript[z, 4] -> 0.955693, Subscript[z, 5] -> -0.254284,
Subscript[z, 6] -> 0.181039, Subscript[z, 7] -> -0.999474,
Subscript[z, 8] -> -0.981127, Subscript[z, 9] -> 0.978473,
Subscript[z, 10] -> 0.999022, Subscript[z, 11] -> -0.0303692,
Subscript[z, 12] -> -0.0391635, Subscript[z, 13] -> -0.151046,
Subscript[z, 14] -> 0.994721, Subscript[z, 15] -> 0.113338,
Subscript[z, 16] -> -0.996332}} *)


Sanity check on those constraints:

Max[Abs[Apply[Subtract, Flatten[constraints], {1}] /. vals]]

(* Out[19]= 7.30056*10^-7 *)


One can do better by setting AccuracyGoal->12 in the FindMinimum[...]. If I do that, the max constraint violation drops to around 10^(-16), minimum staying about the same, at least to all printed digits.

• Off topic but hope you agree MMA is not really updating its optimization capabilities over the last 5-8 years or even more and falling seriously behind as a multipurpose optimization suite. The breakthroughs in algorithms and new class of problems are hardly tractable (e.g using SuperFunctions like FindMinimum ) out of the box. SQP is one such major example. Anything interesting to expect in this direction of making FindMinimum more powerful and compatible to new age demands...? Commented Sep 10, 2013 at 16:05
• Offtopic: Agree. Even some of the third party optimization suites like KNITRO are no longer making Mathematica versions...
– Pam
Commented Sep 10, 2013 at 18:32