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 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}};
share|improve this question
    
Added some sample data to help... –  Pam Sep 9 '13 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. –  Daniel Lichtblau Sep 9 '13 at 22:41
    
Daniel, I have a set of quadratic equations given by data == Transpose [m1.Transpose[θ] + m2.Transpose[θ^2]]; –  Pam Sep 9 '13 at 23:23
    
Edited the original question (f is not re-used any more). –  Pam Sep 9 '13 at 23:26
    
@Daniel: There is not a single objective function but a set of equations subject to a set of constraints... –  Pam Sep 9 '13 at 23:28

2 Answers 2

up vote 3 down vote accepted

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. enter image description here

share|improve this answer
    
Brilliant. This is more than I needed. –  Pam Sep 10 '13 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 Sep 10 '13 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.

share|improve this answer
    
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...? –  PlatoManiac Sep 10 '13 at 16:05
    
Offtopic: Agree. Even some of the third party optimization suites like KNITRO are no longer making Mathematica versions... –  Pam Sep 10 '13 at 18: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.