Unable to compute a rational SOS if the computed solution is almost exact

[PoeticExposition]

Here is an instance:

pvar x y z
p = 4*x^14 - 27*x^12*y^2 + 73*x^10*y^4 + 46*x^10*y^2*z^2 - 35*x^10*z^4 - 104*x^8*y^6 - 77*x^8*y^4*z^2 - 50*x^8*y^2*z^4 + 85*x^8*z^6 + 85*x^6*y^8 + 78*x^6*y^6*z^2 + 7*x^6*y^4*z^4 + 78*x^6*y^2*z^6 - 104*x^6*z^8 - 35*x^4*y^10 - 50*x^4*y^8*z^2 + 7*x^4*y^6*z^4 + 7*x^4*y^4*z^6 - 77*x^4*y^2*z^8 + 73*x^4*z^10 + 46*x^2*y^10*z^2 - 77*x^2*y^8*z^4 + 78*x^2*y^6*z^6 - 50*x^2*y^4*z^8 + 46*x^2*y^2*z^10 - 27*x^2*z^12 + 4*y^14 - 27*y^12*z^2 + 73*y^10*z^4 - 104*y^8*z^6 + 85*y^6*z^8 - 35*y^4*z^10 + 4*z^14;
[Q, ~, ~] = findsos(p, 'rational')

SOSTOOLS says that "Could not compute a rational SOS!", but if we check the output matrix Q without giving the additional input argument 'rational' (i.e., Q = findsos(p);), we'll find that each entry is almost exact indeed. For example,

>> Q(9, 26)
ans =
        4.001
>> Q(24, end)
ans =
        9.999
>> Q(27, 5)
ans =
       -13.99

Well, I have to "rationalize" them manually. Why does the findsos fail to return a symbolic one here?

[djagt]

Hey,

So in establishing a rational SOS decomposition, findsos verifies that the obtained matrix Q is positive semidefinite, by checking whether all the eigenvalues are at least greater than -1e-14. For your problem, numerically computing the eigenvalues of the rounded Q matrix using eig, the smallest eigenvalues turns out to be -3.12e-14, exceeding the implemented threshold. The function findsos takes this as an indicator of the fact that your polynomial does not admit a rational SOS decomposition.

So the problem basically comes down to limited numerical accuracy, as your rational Q matrix is probably positive semidefinite, but we cannot say so for sure by e.g. checking the eigenvalues computed with eig. That being said, in such cases, we should probably set the function to nevertheless return the obtained rational expansion, but display a warning that the smallest (numerically computed) eigenvalue is negative, and thus that the rational Q matrix may not actually be positive semidefinite.

Thanks for pointing this out. I'll push a new version of findsos to GitHub in the near future, informing the user of the small negative eigenvalues. Let me know if you encounter any other issues.

Declan

[PoeticExposition]

So in establishing a rational SOS decomposition, findsos verifies that the obtained matrix Q is positive semidefinite, by checking whether all the eigenvalues are at least greater than -1e-14. For your problem, numerically computing the eigenvalues of the rounded Q matrix using eig, the smallest eigenvalues turns out to be -3.12e-14, exceeding the implemented threshold. The function findsos takes this as an indicator of the fact that your polynomial does not admit a rational SOS decomposition.

Declan

I'm not sure. If I set options.solver = 'sdpa';, the smallest eigenvalues turns out to be 5.7506e-08 (which is greater than -1e-14), but SOSTOOLS still outputs:

Could not compute a rational SOS!
To obtain an exact rational solution, run the function with three output arguments.
Q =
     []

[djagt]

Hey, so of which matrix Q are you checking the eigenvalues? If I run your test with SDPA or another solver, the original (non-rounded) Q matrix indeed has a smallest eigenvalue on the order of +1e-8, but the rounded Q matrix does have a strictly negative eigenvalue. Note that in findsos, the rounded Q matrix is originally stored as Qr.

[djagt]

I'm closing this issue, but feel free to reopen if you feel it hasn't been resolved.

[PoeticExposition]

Sorry for late reply. The above problems are too "large", so it is not easy to display the results here. But yesterday, I found a medium one. The polynomial is:

pvar('a','b','c')
po=@(k)(k*b^2+(a^2+c^2))*(4*b^6+12*(a+c)*b^5+(13*a^2+30*a*c+13*c^2)*b^4+2*(a+3*c)*(3*a+c)*(a+c)*b^3+(a^4+9*a^3*c+12*a^2*c^2+9*a*c^3+c^4)*b^2+a*c*(a+c)*(a-c)^2*b+a^2*c^2*(a-c)^2);%0≤k≤2
op=struct('solver','sdpt3','params',struct('printlevel',0,'scale_data',1,'vers',2),'simplify',1);

SosTools reports that

>> [Q,~,~]=findsos(po(2),'rational',op)
Running simplification process:
Old A size: 169  39
New A size: 169  39

Residual norm: 2.2414e-08

    cpusec: 0.444
      iter: 22
      pinf: 0
      dinf: 0
    numerr: 0
Could not compute a rational SOS!
To obtain an exact rational solution, run the function with three output arguments.
Q =
     []

But if I simply trade b and a, an exact rational solution will be found instead.

>> [Q,~,~]=findsos(subs(po(2),[b;a],[a;b]),'rational',op)
Running simplification process:
Old A size: 169  39
New A size: 169  39

Residual norm: 2.5679e-08

    cpusec: 0.333
      iter: 21
      pinf: 0
      dinf: 0
    numerr: -5
To obtain an exact rational solution, run the function with three output arguments.
Q =
            8           12            4            0           12           16            8            0            4            8            0            0            0
           12           22           12            2           14           24           18            1            0            6            0           -2           -1
            4           12           11            3            0            8           13          1.5           -7           -5            0           -3         -1.5
            0            2            3            1           -2            0            3          0.5           -3           -3            0           -1         -0.5
           12           14            0           -2           22           24            6           -1           12           18            0            2            1
           16           24            8            0           24           40           24           -2            8           24            4            0           -2
            8           18           13            3            6           24           25         -0.5           -5            7            4           -3         -3.5
            0            1          1.5          0.5           -1           -2         -0.5            1         -1.5         -3.5           -1         -0.5            0
            4            0           -7           -3           12            8           -5         -1.5           11           13            0            3          1.5
            8            6           -5           -3           18           24            7         -3.5           13           25            4            3         -0.5
            0            0            0            0            0            4            4           -1            0            4            2            0           -1
            0           -2           -3           -1            2            0           -3         -0.5            3            3            0            1          0.5
            0           -1         -1.5         -0.5            1           -2         -3.5            0          1.5         -0.5           -1          0.5            1

This is strange indeed: findsos(po(2),'rational',op) doesn't work, while findsos(subs(po(2),[b;a],[a;b]),'rational',op) works. Is there any difference? Now I remove the 'rational' flag,

>> [Q,Z,~]=findsos(po(2),op)

Running simplification process:
Old A size: 169  39
New A size: 169  39

Residual norm: 2.2414e-08

    cpusec: 0.3
      iter: 22
      pinf: 0
      dinf: 0
    numerr: 0
Q =
            1            3            2  -2.8207e-06          0.5            3   3.4424e-05           -2  -3.7378e-06           -3           -3         -0.5     -0.99998
            3           11           12            4          1.5           13       8.0001    7.969e-05   -2.694e-05      -4.9999      -6.9999         -1.5           -3
            2           12           22           12            1           18           24           14  -4.9326e-05            6   7.9693e-05     -0.99997           -2
  -2.8207e-06            4           12            8   1.1817e-05       7.9999           16           12  -2.3697e-05       7.9999            4   1.1883e-05  -2.8185e-06
          0.5          1.5            1   1.1817e-05            1         -0.5           -2     -0.99997           -1         -3.5         -1.5    6.705e-10         -0.5
            3           13           18       7.9999         -0.5           25           24            6            4            7      -4.9999         -3.5           -3
   3.4424e-05       8.0001           24           16           -2           24           40           24       3.9999           24       8.0001           -2   3.4381e-05
           -2    7.969e-05           14           12     -0.99997            6           24           22  -4.9339e-05           18           12            1            2
  -3.7378e-06   -2.694e-05  -4.9326e-05  -2.3697e-05           -1            4       3.9999  -4.9339e-05            2            4  -2.7002e-05           -1  -3.7631e-06
           -3      -4.9999            6       7.9999         -3.5            7           24           18            4           25           13         -0.5            3
           -3      -6.9999   7.9693e-05            4         -1.5      -4.9999       8.0001           12  -2.7002e-05           13           11          1.5            3
         -0.5         -1.5     -0.99997   1.1883e-05    6.705e-10         -3.5           -2            1           -1         -0.5          1.5            1          0.5
     -0.99998           -3           -2  -2.8185e-06         -0.5           -3   3.4381e-05            2  -3.7631e-06            3            3          0.5            1
Z = 
  [   a^3*b]
  [ a^2*b^2]
  [   a*b^3]
  [     b^4]
  [   a^3*c]
  [ a^2*b*c]
  [ a*b^2*c]
  [   b^3*c]
  [ a^2*c^2]
  [ a*b*c^2]
  [ b^2*c^2]
  [   a*c^3]
  [   b*c^3]

As you can see, each entry is very close to a rational number. Actually, an exact rational matrix exists:

>> Z.'*arrayfun(@(n)str2num(rat(n,1e-3)),Q)*Z-po(2)
ans = 
  0

But I don't know why SosTools fails to rationalize it. Besides, since the following identity holds:

》pvar k
isequaln(po(k),Z.'*[1 3 2 0 1/2 3 -3/2*k+3 -k 0 -3/2*k -3/2*k -1/2 -k/2;3 k+9 3*k+6 2*k 3/2 2*k+9 -k/2+9 0 0 -5/2*k -7/2*k -3/2 -3/2*k;2 3*k+6 9*k+4 6*k 1 6*k+6 9*k+6 7*k 0 3*k 0 -1 -k;0 2*k 6*k 4*k 0 4*k 8*k 6*k 0 4*k 2*k 0 0;1/2 3/2 1 0 1 -1/2 -2 -1 -1 -7/2 -3/2 0 -1/2;3 2*k+9 6*k+6 4*k -1/2 4*k+17 7/2*k+17 3*k 4 -k/2+8 -5/2*k -7/2 -3/2*k;-3/2*k+3 -k/2+9 9*k+6 8*k -2 7/2*k+17 7*k+26 9*k+6 4 7/2*k+17 -k/2+9 -2 -3/2*k+3;-k 0 7*k 6*k -1 3*k 9*k+6 9*k+4 0 6*k+6 3*k+6 1 2;0 0 0 0 -1 4 4 0 2 4 0 -1 0;-3/2*k -5/2*k 3*k 4*k -7/2 -k/2+8 7/2*k+17 6*k+6 4 4*k+17 2*k+9 -1/2 3;-3/2*k -7/2*k 0 2*k -3/2 -5/2*k -k/2+9 3*k+6 0 2*k+9 k+9 3/2 3;-1/2 -3/2 -1 0 0 -7/2 -2 1 -1 -1/2 3/2 1 1/2;-k/2 -3/2*k -k 0 -1/2 -3/2*k -3/2*k+3 2 0 3 3 1/2 1]*Z)

ans =

  logical

   1

I think that SosTools should be able to find an exact rational solution to po(1). Unfortunately, this time neither findsos(po(1),'rational',op) nor findsos(subs(po(1),[b;a],[a;b]),'rational',op) works. What is the cause?

[djagt]

Thanks for pointing this out, this is an interesting question.

Again, what it seems to come down to is the eigenvalues of the rational matrix Q/Den. For the first test, the rationalized version of the matrix produced by Q=findsos(po(2),op) has a smallest numerical eigenvalue of -1.9716e-14. This is more negative than the specified tolerance -1e-14 in the findsos function. Remarkably, the rationalized version of the matrix produced by Q=findsos(subs(po(2),[a;b],[b;a]),op) has a smallest numerical eigenvalue of -9.0846e-15, less negative than the specified threshold. The function interprets this as the rational expansion of po(2) not being SOS, as the Q matrix has a negative eigenvalue, but the rational expansion of subs(po(2),[a;b],[b;a]) being SOS, as the negative eigenvalue of this Q matrix is just assumed to be a numerical inaccuracy. Presumably these two Q matrices should actually have the same analytical eigenvalues, but I suppose the numerical method for computing these eigenvalues introduces some discrepancy.

I imagine the same issue with the eigenvalues being just slightly too negative also occurs for po(1), though when I ran findsos(po(1),'rational') the function did return an actual rational expansion. Let me know if this continues to be a problem.

I pushed a slightly adjusted version of findsos to GitHub, that should be a bit more accepting of negative eigenvalues, though warning the user of the fact that the matrix may not be positive semidefinite. Please let me know if this does not resolve the issue.