6x6 system of multivariate quadratic equations ... non-negative real solutions

What is the best way (what solver?) to effectively find real non-negative solutions of 6 multivariate quadratic equations (6 variables)?

f_i(x) = 0, i = 1,6, where x = [x1,x2,...,x6] and f_i(x) is the quadratic form with known real coeffs.

[Bruno Luong]

This is a chalenging numerical problem.

The wikipedia article states that Newton or optimization lethods are the current state of the art to solve such system.

By Bezout theorem the number of solutions can go up to 2^6 = 64. All descend numerical methods can get trapped to a one of the local minima (up to 64).

For system of 2 you can read to see how people derive the solutions https://arxiv.org/html/2403.08953v1

You'll appreciate the challenge.

[Michal]

OMG... I just realized how difficult this problem really is! I must completely re-analyze and possibly re-formulate my problem. Thanks for the extremely useful insight!

[Walter Roberson]

It is just possible that if you were to express the equations symbolically, that solve() would be able to come up with solutions. You would probably have to iterate step-by-step and break the system down into parts. It would be a bit tedious... but probably do-able.

[Bruno Luong]

Let's consider a very simple case : If all the quadratic forms are the distance square to 6 given points, then you have a trilateration in 6D space, basically like GPS with 3 more dimensions.

[Walter Roberson]

What is meant here by quadratic equations? Is it the full general form with terms whos total degree does not exceed 2?

f1 = a200000*x(1)^2 + a020000*x(2)^2 + a110000*x(1)*x(2) + a002000*x(3)^2 + a011000*x(2)*x(3) +  a101000*x(1)*x(3) + a000200*x(4)^2 + a001100*x(3)*x(4) + a010100*x(2)*x(4) + a100100*x(1)*x(4) + a000020*x(5)^2 + a000110*x(4)*x(5) + a001010*x(3)*x(5) + a010010*x(2)*x(5) + a100010*x(1)*x(5) + a000002*x(6)^2 + a000011*x(5)*x(6) + a000101*x(4)*x(6) + a0010001*x(3)*x(6) + a0100001*x(2)*x(6) + a100001*x(1)*x(6) + a100000*x(1) + a010000*x(2) + a001000*x(3) + a000100*x(4) + a000010*x(5) + a000001*x(6);

? Or is it something simpler?

[Bruno Luong]

An homogeneuous quadratic equation is of the form

x' * A * x = 0

where the unknown x is n x 1 and A is n x n and given. Without lost of the generality you can assume A is symmetric (hermitian for complex field) .

An non-homogeneuous quadratic equation is of the form

x' * A * x  + b'*x + c = 0

where b is known n x 1 and c is known scalar.

A non-homogeneuous quadratic equation can be considered as homogeneuous quadratic equation in dimension n+1 by the homogenisation as with projective geometry theory. So actually you can considered both forms as your convenience, it just consider the problem of system of quadratic equations in R^5 or R^6 or R^7. I have impression it does simplify a bit the problem.

[Michal]

In my particular problem, c = 0 and the matrix A always contains only one non-zero diagonal element on the ith row/column for f_i(x), i=1,2,...,6 quadratic form. I don't know if these specific properties of quadratic forms in my case can bring any significant simplifications that could lead to a more stable and efficient numerical solution. But in addition, especially in my particular problem, there is an extreme sensitivity to rounding errors. So this is a really challenging numerical problem.

[Walter Roberson]

To confirm, you have the form:

syms x [6 1]
syms A [6 6]
A = diag(diag(A))
x' * A * x

[Bruno Luong]

"A always contains only one non-zero diagonal element on the ith row/column "

To me it means

Ai(j,j) == 0 is true for i ~= j and

false for i == j.

Answer 1

[Bruno Luong]

Assuming the equations are

x' * Ai * x + bi'*x + ci = 0 for i = 1,2, ..., N = 6.

Ai are assumed to be symmetric. If not replace with Asi := 1/2*(Ai + Ai').

You could try to iteratively solve linearized problem; in pseudo code:

x = randn(N,1); % or your solution of previous step, slowly changing as stated
L = zeros(N);
r = zeros(N,1);
notconverge = true;
while notconverge
    for i = 1:N
        L(i,:) = (2*x'*Ai + bi');
        r(i) = -(x'*Ai*x + bi'*x + ci);
    end
    dx = L \ r;
    xold = x;
    x = x + dx;
    x = max(x,0); % since we want solution >= 0
    notconverge  = norm(x-xold,p) > tolx && ...
                   norm(r,q) > tolr; % select p, q, tolx and tolr approproately
end

I guess fmincon, fsolve do somesort of Newton-like or linearization internally. But still worth to investogate. Here the linearization is straight forward and fast to compute. Some intermediate vectors in computing L and r are common and can be shared.

[Michal]

Thanks Bruno!!! This pseudo code looks like a good starting point ... for more thorough investigation.

[Michal]

Because in my case is ci = 0, then r(i)= 0 for x = 0 and i = 1,2, ..., 6

So, in a case of all negative components of x, is then dx = L\r = 0 and convergence of the Newton iterations is broken. Am I right?

[Bruno Luong]

x = 0 is ONE solution (up to 2^6 in the worst case). You have Newton actually converges to a local minima of ONE solution. Bravo.

You need to change

x = randn(N,1);

so as it would converge to a desired solution.

Answer 2

[Aquatris]

fsolve seems to work for most cases, if not lsqnonlin is also an option. here is a general explanation

[Michal]

Frankly, I hope there might be some more specialized solver for quadratic forms.

[John D'Errico]

@Michal

No. There is not. We all want for things that are not possible. Probably more likely to hope for peace in the world. Yeah, right.

Answer 3

[Torsten]

I think there is no specialized solver for a system of quadratic equations. Thus a general nonlinear solver ("fsolve","lsqnonlin","fmincon") is the only chance you have.

[Michal]

The problem is, that I need to solve 6x6 quadratic system many, many ... times (all coefficients are slowly changing) in real-time regime, so speed is the crucial requiremnt.

Any idea how to use these general solvers ("fsolve","lsqnonlin","fmincon") in this regime effectively?

[Torsten]

You shouldn't think about speed at the moment. You are lucky if you get your system solved and if the solution is as expected. A system of quadratic equations is a challenge.

If this works satisfactory, you can save time if you set the solution of the last call to the nonlinear solver to the initial guess of the next call (since you say that your coefficients vary slowly).

[Michal]

Yes... you are definitely right :)