Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Quadratic Cost Function x^T Q x

Subject: Quadratic Cost Function x^T Q x

From: Jason

Date: 20 May, 2010 15:50:21

Message: 1 of 39

Hello,

I have a (quadratic) non-linear problem of the type x^T Q x.

Specifically the non-convex objective function is given by

J = SUM ||x^T Qi x||^2

for i <= 3.

My goal is to find arg min x, i.e. the values for the vector x.

Q is a symmetric 3 x 3 matrix and x is a vector of length 3. This is why i <= 3 since we have 3 parameters (unknowns) in x.

I am currently using lsqnonlin to find a solution. However I am getting trapped in local minima if I don't use a good starting guess for vector x.

My question is, can I use something more 'robust' than lsqnonlin?
Furthermore is there actually a global minimum in my problem?

What would be the best way to find it?

Many thanks.

Subject: Quadratic Cost Function x^T Q x

From: Matt J

Date: 20 May, 2010 16:39:04

Message: 2 of 39

"Jason" <jf203@ic.ac.uk> wrote in message <ht3lnt$92i$1@fred.mathworks.com>...
> Hello,
>
> I have a (quadratic) non-linear problem of the type x^T Q x.
>
> Specifically the non-convex objective function is given by
>
> J = SUM ||x^T Qi x||^2
>
> for i <= 3.
>
> My goal is to find arg min x, i.e. the values for the vector x.
===============

The sum is over i? If so, it is trivial to see that there is a global minimum at x=0.

Subject: Quadratic Cost Function x^T Q x

From: Jason

Date: 20 May, 2010 18:05:20

Message: 3 of 39

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <ht3oj8$gpr$1@fred.mathworks.com>...
> "Jason" <jf203@ic.ac.uk> wrote in message <ht3lnt$92i$1@fred.mathworks.com>...
> > Hello,
> >
> > I have a (quadratic) non-linear problem of the type x^T Q x.
> >
> > Specifically the non-convex objective function is given by
> >
> > J = SUM ||x^T Qi x||^2
> >
> > for i <= 3.
> >
> > My goal is to find arg min x, i.e. the values for the vector x.
> ===============
>
> The sum is over i? If so, it is trivial to see that there is a global minimum at x=0.

Why is there a global minimum at x=0 ??

Maybe I haven't expressed myself right. Qi is a matrix of 3x3, x = column vector of length 3. The result should be 0 if the correct values for x are chosen.

What I mean with the sum is... x^T (Q_1 + Q_2 + ... Q_i) x

Let me show you the code to avoid confusion:

(Note that I need to feed in Q from another function since it changes, Q in this case is a cell of matrices Q.)


function F = opt_routine(x,sns,Q);
for i = 1:sns
    F(i) = (x.')*Q{i}*x;
end
end

We want to find the values for x so:
variable1 = @(x)opt_routine(x,sns,Q);

[x,res] = lsqnonlin(variable1,[] ....)

Ok ?

Subject: Quadratic Cost Function x^T Q x

From: Matt J

Date: 20 May, 2010 18:43:06

Message: 4 of 39

"Jason" <jf203@ic.ac.uk> wrote in message <ht3tl0$mbi$1@fred.mathworks.com>...

> Why is there a global minimum at x=0 ??
>
> Maybe I haven't expressed myself right. Qi is a matrix of 3x3, x = column vector of length 3. The result should be 0 if the correct values for x are chosen.
>
> What I mean with the sum is... x^T (Q_1 + Q_2 + ... Q_i) x
=====================

What you wrote here is not consistent with the code you show below. If the function you are trying to minimize is

f(x) = x^T (Q_1 + Q_2 + ... Q_i) x

Then with A=(Q_1 + Q_2 + ... Q_i), this is a very basic quadratic minimization

f(x)= x'*A*x

If A is positive definite and x is unconstrained, then this function is minimized trivially at x=0. If A is not positive definite and x is unconstrained, this function is not bounded from below, i.e., no global minimum exists. All of this so far is standard theory.

The code below, on the other hand, is as I originally understood your post. According to the code, you want to minimize the function

f(x) = ( x'*Q_1*x )^2 + ( x'*Q_2*x )^2 + ( x'*Q_3*x )^2

This function satisfies f(0)=0. It is also clear that this is a global minimum because
f(x)>=0 for all x (because f(x) is the sum of positive terms).

==================
> Let me show you the code to avoid confusion:
>
> (Note that I need to feed in Q from another function since it changes, Q in this case is a cell of matrices Q.)
>
>
> function F = opt_routine(x,sns,Q);
> for i = 1:sns
> F(i) = (x.')*Q{i}*x;
> end
> end
>
> We want to find the values for x so:
> variable1 = @(x)opt_routine(x,sns,Q);
>
> [x,res] = lsqnonlin(variable1,[] ....)
>
> Ok ?

Subject: Quadratic Cost Function x^T Q x

From: Roger Stafford

Date: 20 May, 2010 18:44:05

Message: 5 of 39

"Jason" <jf203@ic.ac.uk> wrote in message <ht3tl0$mbi$1@fred.mathworks.com>...
>
> Why is there a global minimum at x=0 ??
>
> Maybe I haven't expressed myself right. Qi is a matrix of 3x3, x = column vector of length 3. The result should be 0 if the correct values for x are chosen.
>
> What I mean with the sum is... x^T (Q_1 + Q_2 + ... Q_i) x
>
> Let me show you the code to avoid confusion:
>
> (Note that I need to feed in Q from another function since it changes, Q in this case is a cell of matrices Q.)
>
>
> function F = opt_routine(x,sns,Q);
> for i = 1:sns
> F(i) = (x.')*Q{i}*x;
> end
> end
>
> We want to find the values for x so:
> variable1 = @(x)opt_routine(x,sns,Q);
>
> [x,res] = lsqnonlin(variable1,[] ....)
>
> Ok ?

  You will have to explain yourself better than you have, Jason.

  As things stand, setting all the three elements of x equal to zero guarantees that x'*Qi*x will be the scalar zero, and therefore each component of the sum in J would be zero. Hence J would be zero and it can't be less than zero, so zero would be the minimum. This will remain true until you place some constraint on acceptable x vectors that rules out having all its elements equal to zero.

  If you really mean what you say in x^T (Q_1 + Q_2 + ... Q_i) x, then the eigenvector technique I gave you would still work as applied to the sum of the Q's. However, if you actually mean J is the sum of the squares of the x'*Qi*x scalars, that is an entirely different matter. Which is it you mean?

Roger Stafford

Subject: Quadratic Cost Function x^T Q x

From: Jason

Date: 20 May, 2010 19:08:04

Message: 6 of 39

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <ht3vtl$o1d$1@fred.mathworks.com>...
> "Jason" <jf203@ic.ac.uk> wrote in message <ht3tl0$mbi$1@fred.mathworks.com>...
> >
> > Why is there a global minimum at x=0 ??
> >
> > Maybe I haven't expressed myself right. Qi is a matrix of 3x3, x = column vector of length 3. The result should be 0 if the correct values for x are chosen.
> >
> > What I mean with the sum is... x^T (Q_1 + Q_2 + ... Q_i) x
> >
> > Let me show you the code to avoid confusion:
> >
> > (Note that I need to feed in Q from another function since it changes, Q in this case is a cell of matrices Q.)
> >
> >
> > function F = opt_routine(x,sns,Q);
> > for i = 1:sns
> > F(i) = (x.')*Q{i}*x;
> > end
> > end
> >
> > We want to find the values for x so:
> > variable1 = @(x)opt_routine(x,sns,Q);
> >
> > [x,res] = lsqnonlin(variable1,[] ....)
> >
> > Ok ?
>
> You will have to explain yourself better than you have, Jason.
>
> As things stand, setting all the three elements of x equal to zero guarantees that x'*Qi*x will be the scalar zero, and therefore each component of the sum in J would be zero. Hence J would be zero and it can't be less than zero, so zero would be the minimum. This will remain true until you place some constraint on acceptable x vectors that rules out having all its elements equal to zero.
>
> If you really mean what you say in x^T (Q_1 + Q_2 + ... Q_i) x, then the eigenvector technique I gave you would still work as applied to the sum of the Q's. However, if you actually mean J is the sum of the squares of the x'*Qi*x scalars, that is an entirely different matter. Which is it you mean?
>
> Roger Stafford

Dear Roger,

Thanks for your reply.

I think it was quite clear from the description, and please see the code. lsqnonlin does compute the sum of the squares.

This is the reason why I wrote J = sum || x^T Qi x ||^2

I am referring indeed to the last case you outlined.
Do you have any suggestions for that?

Many thanks

Subject: Quadratic Cost Function x^T Q x

From: Matt J

Date: 20 May, 2010 19:25:21

Message: 7 of 39

"Jason" <jf203@ic.ac.uk> wrote in message <ht41ak$at$1@fred.mathworks.com>...

> I think it was quite clear from the description, and please see the code. lsqnonlin does compute the sum of the squares.
>
> This is the reason why I wrote J = sum || x^T Qi x ||^2
>
> I am referring indeed to the last case you outlined.
> Do you have any suggestions for that?
============

As both Roger and I have suggested already, this function is globally minimized by x=0, so long as x is unconstrained. This is because, by direct evaluation, J=0 for x=0. Also, it is clear that there is no x for which J<0 because J is the sum of positive terms for all x.

It appears as though some constraints are needed to make this problem non-trivial.

Subject: Quadratic Cost Function x^T Q x

From: Roger Stafford

Date: 20 May, 2010 19:33:04

Message: 8 of 39

"Jason" <jf203@ic.ac.uk> wrote in message <ht41ak$at$1@fred.mathworks.com>...
> I think it was quite clear from the description, ......

  No, I disagree with you; your explanation is not at all clear, Jason! You have yet to address the problem of placing constraints on x. Three times that question has been posed to you. Why are you so reticent about that matter? I can make no more suggestions until receiving some statement on your part as to this question.

Roger Stafford

Subject: Quadratic Cost Function x^T Q x

From: Jason

Date: 20 May, 2010 21:22:04

Message: 9 of 39

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <ht42pg$8he$1@fred.mathworks.com>...
> "Jason" <jf203@ic.ac.uk> wrote in message <ht41ak$at$1@fred.mathworks.com>...
> > I think it was quite clear from the description, ......
>
> No, I disagree with you; your explanation is not at all clear, Jason! You have yet to address the problem of placing constraints on x. Three times that question has been posed to you. Why are you so reticent about that matter? I can make no more suggestions until receiving some statement on your part as to this question.
>
> Roger Stafford

Dear Roger,

Without placing constraints on x lsqnonlin runs, and *might* give you the correct result, i.e. the global minimum.

One way to address the problem was to place a constraint on x such that x(1) and x(2) lie on a unit circle. There is no constraint on x(3).

This would reduce the problem to two variables.

if x(1) = cos(alpha) and x(2) = sin(alpha) then the vector x would be:

x = [alpha x(3)].

This works fine using lsqnonlin, however it still does not guarantee that I won't get trapped in a local (suboptimal) minimum.

Think of the vector x as a line in geometrical 2-D space. I want to estimate the parameters of this line, i.e. x(1)*x + x(2)*y + x(3).

Pardon my ignorance, but even without placing any constraints on x you get a non-trivial solution (i.e. not x = 0).

This is why I reformulated my statement, I am not reticent about the matter.

By the way, I have tried fmincon using the nonlinear constraint x(1)^2 + x(2)^2 -1 = 0 but (barring any programming error from my side) did not provide correct results.


My main question is if I can use something different than lsqnonlin which takes a very long time to run and is prone to fall into suboptimal minima.

Best regards,
Jason

Subject: Quadratic Cost Function x^T Q x

From: Matt J

Date: 20 May, 2010 21:47:14

Message: 10 of 39

"Jason" <jf203@ic.ac.uk> wrote in message <ht495s$eq2$1@fred.mathworks.com>...
> "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message

> Think of the vector x as a line in geometrical 2-D space. I want to estimate the parameters of this line, i.e. x(1)*x + x(2)*y + x(3).
>
> Pardon my ignorance, but even without placing any constraints on x you get a non-trivial solution (i.e. not x = 0).
===============

But the cost function for this estimation problem bears no resemblance to the one you posed.


>
> My main question is if I can use something different than lsqnonlin which takes a very long time to run and is prone to fall into suboptimal minima.
=================

You might try running lsqnonlin with the Algorithm option to 'levenberg-marquardt'. My intuition about the default trust-region method is that, even if you initialize in the correct capture basin, the algorithm can crawl out of it into something less optimal. Conversely, Levenberg-Marquardt is, according to my intuition, may be more basin-preserving. Also, if you supply an analytical gradient then, according to doc lsqnonlin, it could be more computationally cheap per iteration.

All of this assumes that you have at least a reasonable guess of the initial solution, but there's no escaping that when it comes to non-convex minimization.

Subject: Quadratic Cost Function x^T Q x

From: Jason

Date: 20 May, 2010 22:07:05

Message: 11 of 39

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <ht4al2$jcc$1@fred.mathworks.com>...
> "Jason" <jf203@ic.ac.uk> wrote in message <ht495s$eq2$1@fred.mathworks.com>...
> > "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message
>
> > Think of the vector x as a line in geometrical 2-D space. I want to estimate the parameters of this line, i.e. x(1)*x + x(2)*y + x(3).
> >
> > Pardon my ignorance, but even without placing any constraints on x you get a non-trivial solution (i.e. not x = 0).
> ===============
>
> But the cost function for this estimation problem bears no resemblance to the one you posed.
>
>
> >
> > My main question is if I can use something different than lsqnonlin which takes a very long time to run and is prone to fall into suboptimal minima.
> =================
>
> You might try running lsqnonlin with the Algorithm option to 'levenberg-marquardt'. My intuition about the default trust-region method is that, even if you initialize in the correct capture basin, the algorithm can crawl out of it into something less optimal. Conversely, Levenberg-Marquardt is, according to my intuition, may be more basin-preserving. Also, if you supply an analytical gradient then, according to doc lsqnonlin, it could be more computationally cheap per iteration.
>
> All of this assumes that you have at least a reasonable guess of the initial solution, but there's no escaping that when it comes to non-convex minimization.

I will try the levenberg-marquardt option, thanks!

As for your inquiry, what do you mean?

The cost function is still the same.

As I said J = min x ( x' Q x).

I am following a principle from computer vision where Q is the matrix representation of a conic (an ellipse). If you set Q = adj(Q) (where adj means adjoint) then x' Q x = 0 means that if this equation is satisfied, then the line parametrized by x is tangent to the ellipse. The additional constraint i have gaven, namely that x(1) and x(2) lie on a unit circle only helps but is not necessary.

Regards,
Jason

Subject: Quadratic Cost Function x^T Q x

From: Matt J

Date: 20 May, 2010 22:56:04

Message: 12 of 39

"Jason" <jf203@ic.ac.uk> wrote in message <ht4bq9$4t9$1@fred.mathworks.com>...

>
> I will try the levenberg-marquardt option, thanks!
>
> As for your inquiry, what do you mean?
>
> The cost function is still the same.
>
> As I said J = min x ( x' Q x).
>
> I am following a principle from computer vision where Q is the matrix representation of a conic (an ellipse). If you set Q = adj(Q) (where adj means adjoint) then x' Q x = 0 means that if this equation is satisfied, then the line parametrized by x is tangent to the ellipse. The additional constraint i have gaven, namely that x(1) and x(2) lie on a unit circle only helps but is not necessary.
=================

You mean x is a vector of homogeneous coordinates!?!

This seems like a really dubious idea. If all this is really just to fit a line, why don't you just use polyfit? It would be much faster and more robust.

At any rate, if you don't constrain homogeneous coordinates in some way, your minimization problem is virtually guaranteed to be ill-conditioned. In the case of the line

x(1)*X+x(2)*Y+x(3)

for any solution x, another solution is c*x for any scalar c. This tends to create problems for optimization algorithms.

Subject: Quadratic Cost Function x^T Q x

From: Jason

Date: 20 May, 2010 23:09:04

Message: 13 of 39

> You mean x is a vector of homogeneous coordinates!?!
>
> This seems like a really dubious idea. If all this is really just to fit a line, why don't you just use polyfit? It would be much faster and more robust.
>
> At any rate, if you don't constrain homogeneous coordinates in some way, your minimization problem is virtually guaranteed to be ill-conditioned. In the case of the line
>
> x(1)*X+x(2)*Y+x(3)
>
> for any solution x, another solution is c*x for any scalar c. This tends to create problems for optimization algorithms.

Yes, x is a vector of homogeneous coordinates!

I agree about what you say: "for any solution x, another solution is c*x for any scalar c." I experienced this, but what can I do about it?

I am reading about how to do this using polyfit, I have no idea how to start, since as I said I want to compute the common tangent to a set of ellipses which is solved if x'Qx = 0.

Say I have a set of ellipses in either matrix form, or of course I can write it out as a polynomial, how do I know if the line I am trying to fit is tangent to these ellipses?

Subject: Quadratic Cost Function x^T Q x

From: Jason

Date: 20 May, 2010 23:30:24

Message: 14 of 39

I have given this some thought, I don't think polyfit is applicable for my purposes.

I want an analytical solution first of all, and secondly I don't know how to compute the common tangent to say 4 ellipses using a method other than the one I have outlined.

Many thanks though.

Subject: Quadratic Cost Function x^T Q x

From: Matt J

Date: 21 May, 2010 03:32:04

Message: 15 of 39

"Jason" <jf203@ic.ac.uk> wrote in message <ht4feg$q4b$1@fred.mathworks.com>...
> > You mean x is a vector of homogeneous coordinates!?!

> I agree about what you say: "for any solution x, another solution is c*x for any scalar c." I experienced this, but what can I do about it?
>
======================

Well, you've already answered that question. By applying the constraint
x(1)^2+x(2)^2=1 for example, not only do you get rid of the continuous 1D space of non-unique solutions c*x, but also the degenerate solution x=0 is excluded from the search.

I'm emphasizing this because you still seem to think the constraint is unnecessary, when in fact it has very important consequences for the performance of the minimization algorithm.


> I am reading about how to do this using polyfit, I have no idea how to start, since as I said I want to compute the common tangent to a set of ellipses which is solved if x'Qx = 0.
===============

No, this is the first time you've explicitly said anything about multiple ellipses. However, now that I understand the problem, I have a few more options to recommend.

First, you can divide the problem into a minimization over two separate regions

(1) The region where we constrain x(3)=0, i.e. where the tangent line passes through the origin.

(2) The region where we constrain x(3)=1, i.e., the region where the tangent line does not pass through the origin. Note that this is equivalent to minimizing over the entire region x(3)~=0 because in homogeneous coordinates x(3)~=0 and x(3)=1 form an equivalence class.

Once you've found the minimum in each case, you can pick the one with the least value of J.

Once we've constrained x(3) either to 0 or 1, there are a few consequences

(I) The cost function J(x) reduces to a 2D function of x(1) and x(2)

(II) It becomes relatively easy to search the 2D space for the global solution by parametrizing the space as
x(1) = t*cosd(alpha)
x(2) = t*sind(alpha)

For each fixed alpha, the cost function J reduces to a 4-th order 1D polynomial in t, which can therefore be minimized rapidly, analytically, and globally in t. I can do this in a loop over a finely sampled range of angles alpha, say alpha=0:0.1:180, until I find a good approximation of the global minimum over all alpha and t and hence over all
x(1) and x(2). This should be fairly quick, considering how easy it is to minimize 1D polynomials

A related approach, but slightly less robust, would be to code your own version of the coordinate descent algorithm. In this algorithm, you alternatingly minimize over each
x(i) while the other x(j) are fixed. Here again J reduces to a quartic polynomial in x(i) which can be globally/analytically minimized in each iteration. Because you are globally minimizing at each step, you have a better than average chance of not getting stuck in a local min. This approach can be applied either to the original 3D problem or one of the above reductions to 2D, although the latter would be better for reasons already discussed.

Subject: Quadratic Cost Function x^T Q x

From: Jason

Date: 21 May, 2010 05:48:04

Message: 16 of 39

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <ht4urj$p85$1@fred.mathworks.com>...
> "Jason" <jf203@ic.ac.uk> wrote in message <ht4feg$q4b$1@fred.mathworks.com>...
> > > You mean x is a vector of homogeneous coordinates!?!
>
> > I agree about what you say: "for any solution x, another solution is c*x for any scalar c." I experienced this, but what can I do about it?
> >
> ======================
>
> Well, you've already answered that question. By applying the constraint
> x(1)^2+x(2)^2=1 for example, not only do you get rid of the continuous 1D space of non-unique solutions c*x, but also the degenerate solution x=0 is excluded from the search.
>
> I'm emphasizing this because you still seem to think the constraint is unnecessary, when in fact it has very important consequences for the performance of the minimization algorithm.
>
>
> > I am reading about how to do this using polyfit, I have no idea how to start, since as I said I want to compute the common tangent to a set of ellipses which is solved if x'Qx = 0.
> ===============
>
> No, this is the first time you've explicitly said anything about multiple ellipses. However, now that I understand the problem, I have a few more options to recommend.
>
> First, you can divide the problem into a minimization over two separate regions
>
> (1) The region where we constrain x(3)=0, i.e. where the tangent line passes through the origin.
>
> (2) The region where we constrain x(3)=1, i.e., the region where the tangent line does not pass through the origin. Note that this is equivalent to minimizing over the entire region x(3)~=0 because in homogeneous coordinates x(3)~=0 and x(3)=1 form an equivalence class.
>
> Once you've found the minimum in each case, you can pick the one with the least value of J.
>
> Once we've constrained x(3) either to 0 or 1, there are a few consequences
>
> (I) The cost function J(x) reduces to a 2D function of x(1) and x(2)
>
> (II) It becomes relatively easy to search the 2D space for the global solution by parametrizing the space as
> x(1) = t*cosd(alpha)
> x(2) = t*sind(alpha)
>
> For each fixed alpha, the cost function J reduces to a 4-th order 1D polynomial in t, which can therefore be minimized rapidly, analytically, and globally in t. I can do this in a loop over a finely sampled range of angles alpha, say alpha=0:0.1:180, until I find a good approximation of the global minimum over all alpha and t and hence over all
> x(1) and x(2). This should be fairly quick, considering how easy it is to minimize 1D polynomials
>
> A related approach, but slightly less robust, would be to code your own version of the coordinate descent algorithm. In this algorithm, you alternatingly minimize over each
> x(i) while the other x(j) are fixed. Here again J reduces to a quartic polynomial in x(i) which can be globally/analytically minimized in each iteration. Because you are globally minimizing at each step, you have a better than average chance of not getting stuck in a local min. This approach can be applied either to the original 3D problem or one of the above reductions to 2D, although the latter would be better for reasons already discussed.

Dear Matt,

Your reply was excellent. This is exactly what I have implemented now.

Many thanks.

Kind regards,
Jason

Subject: Quadratic Cost Function x^T Q x

From: Bruno Luong

Date: 21 May, 2010 05:55:22

Message: 17 of 39

"Jason" <jf203@ic.ac.uk> wrote in message <ht4bq9$4t9$1@fred.mathworks.com>...
>
>
> The cost function is still the same.
>
> As I said J = min x ( x' Q x).
>
> I am following a principle from computer vision where Q is the matrix representation of a conic (an ellipse). If you set Q = adj(Q) (where adj means adjoint) then x' Q x = 0 means that if this equation is satisfied, then the line parametrized by x is tangent to the ellipse. The additional constraint i have gaven, namely that x(1) and x(2) lie on a unit circle only helps but is not necessary.

Can you parametrize the "lines" with the constraints x(1)^2+x^2+x(3)^3 = 1 Rather than just the first two? Seem like it make more symmetry to me.

In any case, such Quadratic programming with norm constraint can be solved by the similar approach I used in my FEX submission:

http://www.mathworks.com/matlabcentral/fileexchange/27596-least-square-with-2-norm-constraint

I believe it is possible to extend the method of norm constraint to semi-norm constraint.

Bruno

Subject: Quadratic Cost Function x^T Q x

From: Matt J

Date: 21 May, 2010 10:33:03

Message: 18 of 39

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ht578a$n72$1@fred.mathworks.com>...

> Can you parametrize the "lines" with the constraints x(1)^2+x^2+x(3)^3 = 1 Rather than just the first two? Seem like it make more symmetry to me.
>
> In any case, such Quadratic programming with norm constraint can be solved by the similar approach I used in my FEX submission:
>
> http://www.mathworks.com/matlabcentral/fileexchange/27596-least-square-with-2-norm-constraint
>
> I believe it is possible to extend the method of norm constraint to semi-norm constraint.
=========

But it is not a quadratic program, Bruno. The objective function is quartic.

Subject: Quadratic Cost Function x^T Q x

From: Bruno Luong

Date: 21 May, 2010 11:09:04

Message: 19 of 39

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <ht5ngv$mot$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message
>
> But it is not a quadratic program, Bruno. The objective function is quartic.

No, in fact they are the same problem: in both cases a quadratic functional is minimized. The formulation in QEP is derived from the Euler Lagrange condition. AFAICS, the calculation *should* be transposed more or less straight forward.

I let OP work on it first, if he get stuck I can give a hand.

Bruno

Subject: Quadratic Cost Function x^T Q x

From: Bruno Luong

Date: 21 May, 2010 11:23:05

Message: 20 of 39

Well, actually I'm confused. The reason I think it's quadratic because I quote (from OP):

'As I said J = min x ( x' Q x).'

OR in another thread

[ If I have a quadratic optimization problem of the type

min(x) x' Q x

Where Q is symmetric but indefinite (3x3) matrix, and x is column vector of length 3.

Fmincon and quadprog could solve this. ]

So what exactly he wants to minimize? It seems Jason confused many people, and even Roger has left the thread.

Bruno

Subject: Quadratic Cost Function x^T Q x

From: Jason

Date: 21 May, 2010 11:59:05

Message: 21 of 39

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ht5qeo$26u$1@fred.mathworks.com>...
> Well, actually I'm confused. The reason I think it's quadratic because I quote (from OP):
>
> 'As I said J = min x ( x' Q x).'
>
> OR in another thread
>
> [ If I have a quadratic optimization problem of the type
>
> min(x) x' Q x
>
> Where Q is symmetric but indefinite (3x3) matrix, and x is column vector of length 3.
>
> Fmincon and quadprog could solve this. ]
>
> So what exactly he wants to minimize? It seems Jason confused many people, and even Roger has left the thread.
>
> Bruno

Bruno, it is quite clearly quadratic.
And it's quite clear that you minimize x, i.e. find the values of x that minimize x' Q x.

If something is not clear in your mind then don't blame me.

Using Matt's approach you can find a solution quickly and efficiently.

He proposed constraining x_3, but you might just consider the subspace x_1 = 1 (which excludes the solution x_1 = 0) so you can combine it with the minimization on a different hyperplane that includes x_1 = 0.

In any case it is possible to obtain a closed form solution and in order to find the global minimum you can find the stationary points of the two cost functions created for each subspace.

Thanks for your help, but I don't think it's in my scope to try it out, I want something from scratch and not hacking a ready code, which I am not even sure is applicable for my purposes.

Best regards,
Jason

Subject: Quadratic Cost Function x^T Q x

From: Matt J

Date: 21 May, 2010 12:03:04

Message: 22 of 39

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ht5qeo$26u$1@fred.mathworks.com>...
 
> So what exactly he wants to minimize? It seems Jason confused many people, and even Roger has left the thread.
==========

This is the objective function

f(x) = ( x'*Q_1*x )^2 + ( x'*Q_2*x )^2 + ( x'*Q_3*x )^2

with some constraint to normalize the magnitude of x, e.g. we have looked at

norm(x)=1
norm(x(1:2))=1
x(3)=1

etc...

Subject: Quadratic Cost Function x^T Q x

From: Bruno Luong

Date: 21 May, 2010 12:24:04

Message: 23 of 39

"Jason" <jf203@ic.ac.uk> wrote in message <ht5si9$ipb$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ht5qeo$26u$1@fred.mathworks.com>...
> > Well, actually I'm confused. The reason I think it's quadratic because I quote (from OP):
> >
> > 'As I said J = min x ( x' Q x).'
> >
> > OR in another thread
> >
> > [ If I have a quadratic optimization problem of the type
> >
> > min(x) x' Q x
> >
> > Where Q is symmetric but indefinite (3x3) matrix, and x is column vector of length 3.
> >
> > Fmincon and quadprog could solve this. ]
> >
> > So what exactly he wants to minimize? It seems Jason confused many people, and even Roger has left the thread.
> >
> > Bruno
>
> Bruno, it is quite clearly quadratic.
> And it's quite clear that you minimize x, i.e. find the values of x that minimize x' Q x.
>
> If something is not clear in your mind then don't blame me.

Well, now Jason, can you tell me the relation between:

1) minimizing x'*Q*x, and
2) x'*Q*x = 0, I quite [ ... common tangent to a set of ellipses which is solved if x'Qx = 0 ]

And earlier on you wrote:

[ ... J = SUM ||x^T Qi x||^2 ... ]

This is clearly quadric functional, not quadratic.

And one more quote, you gave a function handle passing to LSQNONLIN defined by:

[
function F = opt_routine(x,sns,Q);
for i = 1:sns
    F(i) = (x.')*Q{i}*x;
end
end
]

Which LSQNONLIN, will intepret as the cost function as (see http://www.mathworks.com/access/helpdesk/help/toolbox/optim/ug/lsqnonlin.html)

J(x) = sum f(i)^2
= sum ((x.')*Q{i}*x)^2

Which is QUADRIC and not QUADRATIC.

But you state in the last post [ ... it is quite clearly quadratic... ]

There is clearly a confusion in various posts between QUADRIC and QUADRATIC. Who should I blame Jason?

>
> Thanks for your help, but I don't think it's in my scope to try it out, I want something from scratch and not hacking a ready code, which I am not even sure is applicable for my purposes.

If it's quadratic, it's applicable, and FYI it can find ALL the solutions in one shot.

Bruno

Subject: Quadratic Cost Function x^T Q x

From: Bruno Luong

Date: 21 May, 2010 12:28:06

Message: 24 of 39

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <ht5spo$4oc$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ht5qeo$26u$1@fred.mathworks.com>...
>
> > So what exactly he wants to minimize? It seems Jason confused many people, and even Roger has left the thread.
> ==========
>
> This is the objective function
>
> f(x) = ( x'*Q_1*x )^2 + ( x'*Q_2*x )^2 + ( x'*Q_3*x )^2
>

But Jason just tells us [... it is quite clearly quadratic.]!!! And how on earth he can use QUADPROG for this function?

I can't still figure out what Jason really wanted.

Bruno

Subject: Quadratic Cost Function x^T Q x

From: Bruno Luong

Date: 21 May, 2010 12:40:20

Message: 25 of 39

Sorry for the typo, *QUARTIC* and not QUADRIC

Bruno

Subject: Quadratic Cost Function x^T Q x

From: Matt J

Date: 21 May, 2010 13:34:21

Message: 26 of 39

"Jason" <jf203@ic.ac.uk> wrote in message <ht5si9$ipb$1@fred.mathworks.com>...

>
> Bruno, it is quite clearly quadratic.
> And it's quite clear that you minimize x, i.e. find the values of x that minimize x' Q x.
>
> If something is not clear in your mind then don't blame me.
=================

No it is not quadratic. You are trying to find the simultaneous *root* of 3 quadratic functions (one for each ellipse), not their minimum. To do this, you wrote down a cost function J in which you square all the quadratics for each ellipse and add them up. This leads to a 4th order function polynomial in x, not a 2nd order one

Subject: Quadratic Cost Function x^T Q x

From: Jason

Date: 21 May, 2010 15:28:06

Message: 27 of 39

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <ht624t$p6b$1@fred.mathworks.com>...
> "Jason" <jf203@ic.ac.uk> wrote in message <ht5si9$ipb$1@fred.mathworks.com>...
>
> >
> > Bruno, it is quite clearly quadratic.
> > And it's quite clear that you minimize x, i.e. find the values of x that minimize x' Q x.
> >
> > If something is not clear in your mind then don't blame me.
> =================
>
> No it is not quadratic. You are trying to find the simultaneous *root* of 3 quadratic functions (one for each ellipse), not their minimum. To do this, you wrote down a cost function J in which you square all the quadratics for each ellipse and add them up. This leads to a 4th order function polynomial in x, not a 2nd order one

Well it depends on how you define your cost function. lsqnonlin will automatically sum of squares. other methods might not.

In any case, Bruno, I wanted to ask you something else.

I did what you said, i.e. setting x_3 = 1 and x_3 = 0 and then compute two different cost functions.

Also you said we can set x_1 = t * cosd (alpha) and x_2 = t * sind (alpha).

This means that x_1 = x_2 cosd(alpha)/sind(alpha) right ? So if we set beta = cosd(alpha)/sind(alpha) then we have in the first case (x_3 = 0):

J_1 (x_2) = [beta*x_2 x_2 0]*Q*[[beta*x_2 x_2 0]^T

If I do the multiplications correctly then that gives me a polynomial of degree 2 ?! Of course lsqnonlin will sum the squares which would make this degree 4, but what I am saying is correct right?

Also you said I should compute J_1 for each value of alpha from 0:0.1:180 which would make the problem about 1800 iterations for finding x_2. I can then simply find x_1 by computing beta*x_2 and we know that x_3 in this case is 0. Repeat all this for the case that x_3 is 1. Which means 3600 iterations in total.

Have I understood you correctly?

The part which troubles me a bit is when you say: " I can do this in a loop over a finely sampled range of angles alpha, say alpha=0:0.1:180, until I find a good approximation of the global minimum over all alpha and t and hence over all
x(1) and x(2)."

Over all alpha, ok, but over all t ? Haven't I taken t out of the picture here?

Kind regards,
Jason

Subject: Quadratic Cost Function x^T Q x

From: Matt J

Date: 22 May, 2010 15:41:04

Message: 28 of 39

"Jason" <jf203@ic.ac.uk> wrote in message <ht68q6$mcu$1@fred.mathworks.com>...

> Well it depends on how you define your cost function. lsqnonlin will automatically sum of squares. other methods might not.
=====================

That's true, but the cost function you've defined here, and the only one we've been talking about in this thread, is the one used by lsqnonlin. This is for good reason. It would be highly non-standard and take considerable ingenuity to find a simpler cost function for this root-finding problem than what lsqnonlin uses.

 
> In any case, Bruno, I wanted to ask you something else.
==============

Not Bruno. Me.


> I did what you said, i.e. setting x_3 = 1 and x_3 = 0 and then compute two different cost functions.
>
> Also you said we can set x_1 = t * cosd (alpha) and x_2 = t * sind (alpha).
>
> This means that x_1 = x_2 cosd(alpha)/sind(alpha) right ?
================

This manipulation requires a division operation that is only valid for t~=0 and alpha~=0.

It was not my intention that you do this. My intention was that you substitute both
x_1 = t * cosd (alpha) and x_2 = t * sind (alpha)
into the cost function so that, with x_3 and alpha held fixed, the cost function would become a function of t, not x_2.

The result is virtually the same - you end up with a 1D 4th order polynomial to minimize.
However, in your approach, you have to worry about beta being infinite...


 So if we set beta = cosd(alpha)/sind(alpha) then we have in the first case (x_3 = 0):
>
> J_1 (x_2) = [beta*x_2 x_2 0]*Q*[[beta*x_2 x_2 0]^T
>
> If I do the multiplications correctly then that gives me a polynomial of degree 2 ?! Of course lsqnonlin will sum the squares which would make this degree 4, but what I am saying is correct right?
===================

You shouldn't be using lsqnonlin for this. lsqnonlin does not know how to minimize anything analytically, even a simple 1D polynomial. You should be writing your own routine to do this.

To keep this discussion clear, it would also be a good idea if you stop treating lsqnonlin as the thing that's defining the cost function. WE are defining the cost function for this problem and it is the multi-variable 4th order polynomial

J(x)=sum_i (x' * Q{i} *x)^2

There is no other obvious and more tractable alternative to this cost function for the task you've described.


> Also you said I should compute J_1 for each value of alpha from 0:0.1:180 which would make the problem about 1800 iterations for finding x_2. I can then simply find x_1 by computing beta*x_2 and we know that x_3 in this case is 0. Repeat all this for the case that x_3 is 1. Which means 3600 iterations in total.
===============

One other thing. It has occured to me that you don't really have to do this for the case x_3=0

In this case, the remaining unknowns x_1 and x_2 can also be normalized with impunity so that x_1=1 or x_2=1 as long as one of these is non-zero. So, you can subdivide the case x_3=0 into 2 more subcases

(1) x_1=1, x_3=0
(2) x_2=1, x_3=0

In each of these sub-cases, the cost function reduces to a 1D polynomial again, which you can minimize trivially.

Subject: Quadratic Cost Function x^T Q x

From: Jason

Date: 22 May, 2010 23:42:04

Message: 29 of 39

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <ht8tug$sjj$1@fred.mathworks.com>...
> "Jason" <jf203@ic.ac.uk> wrote in message <ht68q6$mcu$1@fred.mathworks.com>...
>
> > Well it depends on how you define your cost function. lsqnonlin will automatically sum of squares. other methods might not.
> =====================
>
> That's true, but the cost function you've defined here, and the only one we've been talking about in this thread, is the one used by lsqnonlin. This is for good reason. It would be highly non-standard and take considerable ingenuity to find a simpler cost function for this root-finding problem than what lsqnonlin uses.
>
>
> > In any case, Bruno, I wanted to ask you something else.
> ==============
>
> Not Bruno. Me.
>
>
> > I did what you said, i.e. setting x_3 = 1 and x_3 = 0 and then compute two different cost functions.
> >
> > Also you said we can set x_1 = t * cosd (alpha) and x_2 = t * sind (alpha).
> >
> > This means that x_1 = x_2 cosd(alpha)/sind(alpha) right ?
> ================
>
> This manipulation requires a division operation that is only valid for t~=0 and alpha~=0.
>
> It was not my intention that you do this. My intention was that you substitute both
> x_1 = t * cosd (alpha) and x_2 = t * sind (alpha)
> into the cost function so that, with x_3 and alpha held fixed, the cost function would become a function of t, not x_2.
>
> The result is virtually the same - you end up with a 1D 4th order polynomial to minimize.
> However, in your approach, you have to worry about beta being infinite...
>
>
> So if we set beta = cosd(alpha)/sind(alpha) then we have in the first case (x_3 = 0):
> >
> > J_1 (x_2) = [beta*x_2 x_2 0]*Q*[[beta*x_2 x_2 0]^T
> >
> > If I do the multiplications correctly then that gives me a polynomial of degree 2 ?! Of course lsqnonlin will sum the squares which would make this degree 4, but what I am saying is correct right?
> ===================
>
> You shouldn't be using lsqnonlin for this. lsqnonlin does not know how to minimize anything analytically, even a simple 1D polynomial. You should be writing your own routine to do this.
>
> To keep this discussion clear, it would also be a good idea if you stop treating lsqnonlin as the thing that's defining the cost function. WE are defining the cost function for this problem and it is the multi-variable 4th order polynomial
>
> J(x)=sum_i (x' * Q{i} *x)^2
>
> There is no other obvious and more tractable alternative to this cost function for the task you've described.
>
>
> > Also you said I should compute J_1 for each value of alpha from 0:0.1:180 which would make the problem about 1800 iterations for finding x_2. I can then simply find x_1 by computing beta*x_2 and we know that x_3 in this case is 0. Repeat all this for the case that x_3 is 1. Which means 3600 iterations in total.
> ===============
>
> One other thing. It has occured to me that you don't really have to do this for the case x_3=0
>
> In this case, the remaining unknowns x_1 and x_2 can also be normalized with impunity so that x_1=1 or x_2=1 as long as one of these is non-zero. So, you can subdivide the case x_3=0 into 2 more subcases
>
> (1) x_1=1, x_3=0
> (2) x_2=1, x_3=0
>
> In each of these sub-cases, the cost function reduces to a 1D polynomial again, which you can minimize trivially.

Matt, yes, sorry about this!

Re: My intention was that you substitute both
> x_1 = t * cosd (alpha) and x_2 = t * sind (alpha)
> into the cost function

Yes, just realized that before you answered. Got you now!

However I do not understand your last point:

> One other thing. It has occured to me that you don't really have to do this for the case x_3=0
>
> In this case, the remaining unknowns x_1 and x_2 can also be normalized with impunity so that x_1=1 or x_2=1 as long as one of these is non-zero. So, you can subdivide the case x_3=0 into 2 more subcases.

Could you elaborate this a little bit please, I'm not sure I get you. You are basically saying than when x_3 = 0 then x_1 and x_2 ~= 0 ?

I don't really understand.

Thanks!

Subject: Quadratic Cost Function x^T Q x

From: Matt J

Date: 23 May, 2010 14:24:03

Message: 30 of 39

"Jason" <jf203@ic.ac.uk> wrote in message <ht9q4c$jbj$1@fred.mathworks.com>...

> Could you elaborate this a little bit please, I'm not sure I get you. You are basically saying than when x_3 = 0 then x_1 and x_2 ~= 0 ?
>
> I don't really understand.
=========

It's the same equivalence class consideration associated with homogeneous coordinates that we already said made the space of solutions with x3~=0 equivalent to the space of solutions with x3=1.

Once you restrict the search to the space x_3=0, it means you are searching for a solution of the form (x1, x2, 0). Suppose we restrict this search further to the subspace of solutions of the form (x1,x2,0) and for which x1~=0. Each of the solutions in this subspace are equivalent to a solution of the form
(1,x2/x1,0) or more simply (1,z,0). Therefore, you can just as well restrict the search to the subspace (1,x2,0)

You can do the same thing for the subspace {x2~=0, x3=0}
It is equivalent to the subspace {x2=1, x3=0}

Subject: Quadratic Cost Function x^T Q x

From: Jason

Date: 19 Jul, 2010 18:35:05

Message: 31 of 39

I seem to have hit a brick-wall, when I initially thought the problem was solved.

Just to re-iterate, we are facing a system of the type (x^T Q x)^2 for which x = [x1,x2,x3] and Q is a 3x3 matrix representing a conic (or more specifically an ellipse). Note that we are dealing with homogenous coordinates.

The idea was to find the common tangent to a set of ellipses Q1, Q2, Q3, Q4 ... which is done by computing the adjoint of sum{Q1,...,Q4} and then solving for the parameters x1,x2,x3.

We did this (as proposed by Matt) by decomposing the problem into 3 cases:

Case 1: Set x = [x1,x2,1]

Case 2: Set x = [x1,1,0]

Case 3: Set x = [1,x2,0]

Obtaining the solutions for Cases 2 and 3 is trivial as we are dealing with a 1-D polynomial for which we can compute the roots and hence obtain our solution.

Case 1 is a bit more complex since we are dealing with a fourth order polynomial in both x1 and x2.

For this case we have our initial function f = ( [x1,x2,1]^T * Q * [x1,x2,1] )^2.
Note that Q = [a,b,d;b,c,e;d,e,f].

Hence f = (a*x1^2 + 2b*x1*x2 + 2d*x1 + c*x2^2 + 2e*x2 + f)^2

What I do next is compute the partial derivative of f with respect to x1 and the partial derivative of f with respect to x2. I then set the first partial derivative to 0 to compute the roots:

df/dx1 = 0.

I am left with some roots (potentially complex) which I substitute for x1 into the partial derivative of f with respect to x2, i.e. into df/dx2.

I set this to zero as well to get my solutions for x2, i.e. I do

df/dx2 (x1,1,...,x1,N) = 0.

I am left with an expression for x2 which is only dependent on Q, i.e. on known parameters a,b,c,d,e,f !

These solutions are potentially complex as well.

I then take these results for x2 and resubstitute into the original roots i found for df/dx1 = 0 to compute values x1.

This is pretty standard practice for computing the global minima of a function depending on two variables right?


The problem is, my results are wrong. And I don't know why.

Is there something wrong with my thinking? Did I omit something. Is something not true when dealing with potentially complex roots?

I can be more specific if you like!

Many thanks

Subject: Quadratic Cost Function x^T Q x

From: Matt J

Date: 19 Jul, 2010 20:17:05

Message: 32 of 39

"Jason" <jf203@ic.ac.uk> wrote in message <i225sp$gej$1@fred.mathworks.com>...

> The idea was to find the common tangent to a set of ellipses Q1, Q2, Q3, Q4 ... which is done by computing the adjoint of sum{Q1,...,Q4} and then solving for the parameters x1,x2,x3.
==========

No, you cannot approach the problem by combining all the Q{i} into one matrix. The problem as we've posed it is to minimize,

J(x)=sum_i (x' * adjoint(Q{i}) *x)^2

This is NOT equivalent to the function (x' * sum_i (adjoint(Q{i}) *x)^2

If that had been true, it would mean that the problem of finding a line uniquely tangent to multiple ellipses is equivalent to finding a line intersecting the single ellipse sum_i (Q{i}). But since a single ellipse never has a unique tangent line, you know it can't be equivalent to a multi-ellipse case in which the solution is unique.


> What I do next is compute the partial derivative of f with respect to x1 and the partial derivative of f with respect to x2. I then set the first partial derivative to 0 to compute the roots:
>
> df/dx1 = 0.
>
> I am left with some roots (potentially complex)
=================

You wouldn't get numerical values for these roots. At best, you would get VERY complicated expressions for the roots in terms of the unknown x2.




>which I substitute for x1 into the partial derivative of f with respect to x2, i.e. into df/dx2.
>
> I set this to zero as well to get my solutions for x2, i.e. I do
>
> df/dx2 (x1,1,...,x1,N) = 0.
>
> I am left with an expression for x2 which is only dependent on Q, i.e. on known parameters a,b,c,d,e,f !
===================
 
But because the expressions substituted in for x1 are very complicated functions of x2, this equation should be a VERY hard root finding problem. How did you solve it?


> These solutions are potentially complex as well.
>
> I then take these results for x2 and resubstitute into the original roots i found for df/dx1 = 0 to compute values x1.
>
> This is pretty standard practice for computing the global minima of a function depending on two variables right?
===============

If you can simultaneously solve the equations df/dxi=0, then yes, but often those equations are too complicated to solve analytically (as they are here) and need to be solved by more indirect means. That's why the approaches I recommended to you in earlier posts were more elaborate.

Subject: Quadratic Cost Function x^T Q x

From: Matt J

Date: 19 Jul, 2010 20:29:35

Message: 33 of 39

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i22bs1$bsh$1@fred.mathworks.com>...
>The problem as we've posed it is to minimize,
>
> J(x)=sum_i (x' * adjoint(Q{i}) *x)^2
>
> This is NOT equivalent to the function (x' * sum_i (adjoint(Q{i}) *x)^2
=======

Should have been ( x' * adjoint(sum_i Q{i}) * x )^2

Subject: Quadratic Cost Function x^T Q x

From: Jason

Date: 19 Jul, 2010 21:04:06

Message: 34 of 39

No, you cannot approach the problem by combining all the Q{i} into one matrix. The problem as we've posed it is to minimize,

J(x)=sum_i (x' * adjoint(Q{i}) *x)^2

This is NOT equivalent to the function (x' * sum_i (adjoint(Q{i}) *x)^2

==============================================

Point taken!

You wouldn't get numerical values for these roots. At best, you would get VERY complicated expressions for the roots in terms of the unknown x2.

===============================================

Yes, exactly. It will be a complicated expression in terms of the unknown x2. However you can solve it using symlinks and the matlab 'solve' command.

For example (NOTE THAT x1 and x2 are replaced by l1 and l2 in my code!):

J1 = (a*l1^2 + 2*b*l1*l2 + 2*d*l1 + c*l2^2 + 2*e*l2 + f)^2

The first derivative in terms of l1:
dJl1 = 4*(d + a*l1 + b*l2)*(a*l1^2 + 2*b*l1*l2 + 2*d*l1 + c*l2^2 + 2*e*l2 + f)

The roots of this:

root1 = -(d + b*l2)/a
root2 = -(d + b*l2 - (b^2*l2^2 + 2*b*d*l2 + d^2 - a*c*l2^2 - 2*a*e*l2 - a*f)^(1/2))/a
root3 = -(d + b*l2 + (b^2*l2^2 + 2*b*d*l2 + d^2 - a*c*l2^2 - 2*a*e*l2 - a*f)^(1/2))/a

Solutions for l2 based on the first root:

l2_1 = (b*d - a*e + a*((f*b^2 - 2*b*d*e + c*d^2 + a*e^2 - a*c*f)/a)^(1/2))/(a*c - b^2)
l2_2 = -(a*e - b*d)/(a*c - b^2)
l2_3 = -(a*e - b*d + a*((f*b^2 - 2*b*d*e + c*d^2 + a*e^2 - a*c*f)/a)^(1/2))/(a*c - b^2)

Solutions for l1 by substituting back into root1:

l1_1 = -(d + (b*(b*d - a*e + a*((f*b^2 - 2*b*d*e + c*d^2 + a*e^2 - a*c*f)/a)^(1/2)))/(a*c - b^2))/a
l1_2 = -(d - (b*(a*e - b*d))/(a*c - b^2))/a
l1_3 = -(d - (b*(a*e - b*d + a*((f*b^2 - 2*b*d*e + c*d^2 + a*e^2 - a*c*f)/a)^(1/2)))/(a*c - b^2))/a

etc...

Is there anything wrong with that?

How else could I approach the problem then?

And wouldn't it be possible to just correct the first inconsistency you presented by summing over i in all the cases but following the same methodology ?

Regards,
Jason

Subject: Quadratic Cost Function x^T Q x

From: Matt J

Date: 20 Jul, 2010 15:04:05

Message: 35 of 39

"Jason" <jf203@ic.ac.uk> wrote in message <i22ek6$90o$1@fred.mathworks.com>...

> Yes, exactly. It will be a complicated expression in terms of the unknown x2. However you can solve it using symlinks and the matlab 'solve' command.
>
> For example (NOTE THAT x1 and x2 are replaced by l1 and l2 in my code!):
>
> J1 = (a*l1^2 + 2*b*l1*l2 + 2*d*l1 + c*l2^2 + 2*e*l2 + f)^2
>
> The first derivative in terms of l1:
> dJl1 = 4*(d + a*l1 + b*l2)*(a*l1^2 + 2*b*l1*l2 + 2*d*l1 + c*l2^2 + 2*e*l2 + f)
>
> The roots of this:
>
> root1 = -(d + b*l2)/a
> root2 = -(d + b*l2 - (b^2*l2^2 + 2*b*d*l2 + d^2 - a*c*l2^2 - 2*a*e*l2 - a*f)^(1/2))/a
> root3 = -(d + b*l2 + (b^2*l2^2 + 2*b*d*l2 + d^2 - a*c*l2^2 - 2*a*e*l2 - a*f)^(1/2))/a
>
> Solutions for l2 based on the first root:
>
> l2_1 = (b*d - a*e + a*((f*b^2 - 2*b*d*e + c*d^2 + a*e^2 - a*c*f)/a)^(1/2))/(a*c - b^2)
> l2_2 = -(a*e - b*d)/(a*c - b^2)
> l2_3 = -(a*e - b*d + a*((f*b^2 - 2*b*d*e + c*d^2 + a*e^2 - a*c*f)/a)^(1/2))/(a*c - b^2)
>
> Solutions for l1 by substituting back into root1:
>
> l1_1 = -(d + (b*(b*d - a*e + a*((f*b^2 - 2*b*d*e + c*d^2 + a*e^2 - a*c*f)/a)^(1/2)))/(a*c - b^2))/a
> l1_2 = -(d - (b*(a*e - b*d))/(a*c - b^2))/a
> l1_3 = -(d - (b*(a*e - b*d + a*((f*b^2 - 2*b*d*e + c*d^2 + a*e^2 - a*c*f)/a)^(1/2)))/(a*c - b^2))/a
>
> etc...
>
> Is there anything wrong with that?
===========

It is strange to me that for each fixed root1,2,3 the symbolic solver would obtain only a single corresponding l1_1,2,3.

Prior to substituting root1,2,3 into dJl2, you have a cubic polynomial in l2 with potentially multiple real roots for each fixed l1. Are you saying the substitution of root1,2,3 reduced the number of potential roots to one?
 

> How else could I approach the problem then?
=========


We already discussed an alternative in Message #15 and #28


> And wouldn't it be possible to just correct the first inconsistency you presented by summing over i in all the cases but following the same methodology ?
========

Yes, even after the correction, you would obtain 4th order polynomials for J(l1,l2) and 3rd order polynomials for its derivatives. Any method applicable to that scenario should still be applicable after the correction.

Subject: Quadratic Cost Function x^T Q x

From: Jason

Date: 21 Jul, 2010 13:25:05

Message: 36 of 39

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i24dt5$3op$1@fred.mathworks.com>...
> "Jason" <jf203@ic.ac.uk> wrote in message <i22ek6$90o$1@fred.mathworks.com>...
>
> > Yes, exactly. It will be a complicated expression in terms of the unknown x2. However you can solve it using symlinks and the matlab 'solve' command.
> >
> > For example (NOTE THAT x1 and x2 are replaced by l1 and l2 in my code!):
> >
> > J1 = (a*l1^2 + 2*b*l1*l2 + 2*d*l1 + c*l2^2 + 2*e*l2 + f)^2
> >
> > The first derivative in terms of l1:
> > dJl1 = 4*(d + a*l1 + b*l2)*(a*l1^2 + 2*b*l1*l2 + 2*d*l1 + c*l2^2 + 2*e*l2 + f)
> >
> > The roots of this:
> >
> > root1 = -(d + b*l2)/a
> > root2 = -(d + b*l2 - (b^2*l2^2 + 2*b*d*l2 + d^2 - a*c*l2^2 - 2*a*e*l2 - a*f)^(1/2))/a
> > root3 = -(d + b*l2 + (b^2*l2^2 + 2*b*d*l2 + d^2 - a*c*l2^2 - 2*a*e*l2 - a*f)^(1/2))/a
> >
> > Solutions for l2 based on the first root:
> >
> > l2_1 = (b*d - a*e + a*((f*b^2 - 2*b*d*e + c*d^2 + a*e^2 - a*c*f)/a)^(1/2))/(a*c - b^2)
> > l2_2 = -(a*e - b*d)/(a*c - b^2)
> > l2_3 = -(a*e - b*d + a*((f*b^2 - 2*b*d*e + c*d^2 + a*e^2 - a*c*f)/a)^(1/2))/(a*c - b^2)
> >
> > Solutions for l1 by substituting back into root1:
> >
> > l1_1 = -(d + (b*(b*d - a*e + a*((f*b^2 - 2*b*d*e + c*d^2 + a*e^2 - a*c*f)/a)^(1/2)))/(a*c - b^2))/a
> > l1_2 = -(d - (b*(a*e - b*d))/(a*c - b^2))/a
> > l1_3 = -(d - (b*(a*e - b*d + a*((f*b^2 - 2*b*d*e + c*d^2 + a*e^2 - a*c*f)/a)^(1/2)))/(a*c - b^2))/a
> >
> > etc...
> >
> > Is there anything wrong with that?
> ===========
>
> It is strange to me that for each fixed root1,2,3 the symbolic solver would obtain only a single corresponding l1_1,2,3.
>
> Prior to substituting root1,2,3 into dJl2, you have a cubic polynomial in l2 with potentially multiple real roots for each fixed l1. Are you saying the substitution of root1,2,3 reduced the number of potential roots to one?
>

I have only outlined the case of substituting root1 into dJl2 and setting to zero.

In fact if you do the same for root2 and root3 you obtain 9 roots in total (from which 2 can be excluded) meaning you are left with 7 in total, i.e. 7 estimates for l1 and consequently 7 estimates for l2.

Since I wrongly placed the sum inside rather than outside I am missing some terms which I need to update now. But the solution should still be straightforward, and most importantly analytical and closed form.

I know we have discussed about parametrizing l1 and l2 using cos and sin. However this is not an option I'm afraid.

Do you see any problems with the way I am about to proceed ?

Subject: Quadratic Cost Function x^T Q x

From: Matt J

Date: 21 Jul, 2010 13:37:04

Message: 37 of 39

"Jason" <jf203@ic.ac.uk> wrote in message <i26sfh$spl$1@fred.mathworks.com>...

>
> I know we have discussed about parametrizing l1 and l2 using cos and sin. However this is not an option I'm afraid.
>
> Do you see any problems with the way I am about to proceed ?

Not if the symbolic toolbox is able to solve the optimality equations, as you say. The only problem I see is the unclear reasoning why parametrizing l1 and l2 with sines and cosines is "not an option".

Subject: Quadratic Cost Function x^T Q x

From: Jason

Date: 21 Jul, 2010 13:54:05

Message: 38 of 39

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i26t60$f5i$1@fred.mathworks.com>...
> "Jason" <jf203@ic.ac.uk> wrote in message <i26sfh$spl$1@fred.mathworks.com>...
>
> >
> > I know we have discussed about parametrizing l1 and l2 using cos and sin. However this is not an option I'm afraid.
> >
> > Do you see any problems with the way I am about to proceed ?
>
> Not if the symbolic toolbox is able to solve the optimality equations, as you say.

This I am not 100% sure of, but am testing it now. I did come across some problems for which I had to do some workarounds.



Parametrizing l1 and l2 is not an option because I would have to do an exhaustive search over all angles. Furthermore this will only be an estimate, and yes, it can be improved by using a finer grid once a solution has been found.

I am not adverse to this, and might actually give it a try if everything else fails. However I am worried about speed, and also not having an exact solution, when theoretically you can obtain a perfect solution, i.e. analytically obtain the global minimum.

That's why!

Many thanks for your insightful feedback, I will let you know how it went!

Subject: Quadratic Cost Function x^T Q x

From: Matt J

Date: 21 Jul, 2010 14:08:05

Message: 39 of 39

"Jason" <jf203@ic.ac.uk> wrote in message <i26u5t$jk6$1@fred.mathworks.com>...

>
> Parametrizing l1 and l2 is not an option because I would have to do an exhaustive search over all angles. Furthermore this will only be an estimate, and yes, it can be improved by using a finer grid once a solution has been found.
=============

No, that's not the way I was intending you go about it. You would search on a grid of angles that is both coarse enough to make the search very fast search, but fine enough that the search gets you very close to the solution. Depending on the size of your ellipses, the grid size might be 1 or 2 degrees or even more.

Once you have the result of this search, you would use it to initialize lsqnonlin and converge to a more exact result. The problems you were initially having with lsqnonlin would be greatly diminished here because now you have a very good initial guess of the global minimum. That will make the residual work done by lsqnonlin both fast and robust against stuck at local minima.

> I am not adverse to this, and might actually give it a try if everything else fails. However I am worried about speed, and also not having an exact solution, when theoretically you can obtain a perfect solution, i.e. analytically obtain the global minimum.
================

I'm still rather astounded that this can be solved exactly, but yes, an analytical solution is better.

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us