%**************************************************************************
% 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. |
Documentation > NCSOStoolsdemo >