%**************************************************************************
% Rationalization *
%**************************************************************************
% In this section we will look at the module RprojRldlt which tries to find
% an exact rational positive semidefinite solution of the problem A*X=b
% from the floating point solution XX. The algorithm tries to find an
% appropiate rational approximation of the matrix XX which is after the
% projection on the hyperplane A*X=b a positive semidefinite matrix.
% Module has also an option to try to compute an exact rational LDU
% decomposition of a rational positive semidefinite solution of the problem
% A*XX=b. Otherwise positive semidefiniteness relies on numerical
% computation of eigenvalues.
% This is usefull whenever we want to find a rational Gram matrix an
% therefore a rational sum of hermitian squares decomposition, i.e. in
% module NCsos, NCmin, NCcycSos, NCcycMin ...
% -------------------------------------------------------------------------
% First we construct some symbolic noncommuting variables and polynomial
% -------------------------------------------------------------------------
>> NCvars x y
>> f = 1 + 3*x^2*y*x^2*y + 7*x^2*y^2*x^2 + 2*x*y*x + 6*x*y*x*y*x^2;
% -------------------------------------------------------------------------
% With NCcycSos we can get numerical solution and therefore numerical
% (inexact) sum of hermitian squares decomposition. This will be our
% starting point in finding exact solution.
% For the rationalization is in most cases better to take objective
% function for the SDP solver to be zero; so we set this in parameter.
% -------------------------------------------------------------------------
>> params.obj = 0; [IsCycEq,XX,base,sohs,g,SDP_data] = NCcycSos(f, params)
***** NCSOStools: module NCcycSos started *****
Input polynomial has (max) degree 6 and min degree 0.
Detected 5 monomials in 2 variables.
There are 127 monomials in 2 variables of degree at most 6.
Preprocessing the input polynomial ...
Using alpha degree to construct the monomial vector.
Keeping 4 monomials in the monomial vector.
Computing cyclically equivalent products ... done.
Preparing linear constraints ...
Number of linear constraints: 5.
Starting SDP solver ...
...
Residual norm: 4.6902e-014
***** ALERT: SDP solver ran into minor numerical problems. *****
***** Nevertheless, according to feasratio (0.8878) given solution MIGHT be ok. *****
Computing SOHS decomposition ... done.
Found SOHS decomposition with 4 factors.
WARNING! SDP solver returned some numerical problems. Check messages!
IsCycEq = 1
XX = 1.0000 0.0439 0.5087 0.4474
0.0439 2.1604 0.6947 -1.7065
0.5087 0.6947 6.4130 2.3053
0.4474 -1.7065 2.3053 4.8396
base = ''
'x*x*y'
'x*y*x'
'y*x*x'
sohs = 1+0.0439*x^2*y+0.5087*x*y*x+0.4474*y*x^2
1.469*x^2*y+0.4577*x*y*x-1.175*y*x^2
2.438*x*y*x+1.073*y*x^2
1.452*y*x^2
g = 1+0.4913*x^2*y - 1.7065222*x^2*y*x^2*y + 2.3052978*x^2*y*x*y*x + 4.8395461*x^2*y^2*x^2 + 1.0174*x*y*x + 0.69478477*x*y*x^3*y + 6.4130842*x*y*x^2*y*x + 2.3052978*x*y*x*y*x^2 + 0.4913*y*x^2 + 2.1604759*y*x^4*y + 0.69478477*y*x^3*y*x - 1.7065222*y*x^2*y*x^2
SDP_data =
A: [5x16 double]
b: [5x1 double]
C: [16x1 double]
K: [1x1 struct]
pars: [1x1 struct]
% -------------------------------------------------------------------------
% Now we can use the algorithm in RprojRldlt to try to find an appropiate
% rational approximation of the matrix XX which is after the projection on
% the hyperplane A*X=b a positive semidefinite matrix; if we succeed in
% finding rational Gram matrix we can therefore get a rational sum of
% hermitian squares decomposition. The next step depends heavily on the
% 'quality' of the Gram matrix XX returned by the SDP solver. Here we
% measure 'quality' as large as possible minimal eigenvalue of XX.
% -------------------------------------------------------------------------
>> Xrat = RprojRldlt(XX,SDP_data,true); Xrat{1}, Xrat{2}
***** NCSOStools: module RProjRldlt started *****
delta (min eig X): 7.373264e-001
epsilon (residual norm): 4.690190e-014
i : 1/2^i tau_i min_eig
0 : 1.00E+000 8.79E-001 4.49E-001
Rational projection and LDLT decomposition found.
***** NCSOStools: module RProjRldlt completed successfully *****
ans = 1 0 1 0
0 2 1 -5
1 1 19 2
0 -5 2 5
ans = 1 1 1 1
1 1 1 3
1 1 3 1
1 3 1 1
% -------------------------------------------------------------------------
% Now we got an exact rational Gram matrix Xrat where Xrat{1} are
% numerators and Xrat{2} adherent denominators. If we also want a rational
% sum of hermitian squares decomposition we need to slightly change a call
% to the module RprojRldlt to get a rational LDU decomposition of Xrat,
% more precise we get a rational decomposition of the form
% Xrat=(P*L)*D*(P*L)^T.
% -------------------------------------------------------------------------
>> [Xrat,L,D,P,err] = RprojRldlt(XX,SDP_data,true);
***** NCSOStools: module RProjRldlt started *****
delta (min eig X): 7.373264e-001
epsilon (residual norm): 4.690190e-014
i : 1/2^i tau_i min_eig
0 : 1.00E+000 8.79E-001 4.49E-001
Rational projection and LDLT decomposition found.
***** NCSOStools: module RProjRldlt completed successfully *****
>> base
base = ''
'x*x*y'
'x*y*x'
'y*x*x'
>> P
P = 0 0 0 1
0 0 1 0
1 0 0 0
0 1 0 0
>> D{1}, D{2}
ans = 19 0 0 0
0 83 0 0
0 0 704 0
0 0 0 509
ans = 3 1 1 1
1 19 1 1
1 1 747 1
1 1 1 704
>> L{1}, L{2}
ans = 1 0 0 0
6 1 0 0
3 -113 1 0
3 -6 -225 1
ans = 1 1 1 1
19 1 1 1
19 249 1 1
19 83 704 1
% If RprojRldlt was succesfull we gain an exact rational in the following
% way:
% - first we permute a base vector according to the permutation matrix P
% - each column of the matrix L corresponds to the coefficients of the one
% hermitian square in the decomposition multiplied by the corresponding
% element in the matrix D.