%**************************************************************************
% Sum of hermitian squares (SOHS) related                                 *
%**************************************************************************

% -------------------------------------------------------------------------
% First we construct some symbolic noncommuting variables and some
% polynomials
% -------------------------------------------------------------------------

>> NCvars x y z
>> f = y*x^2*y - y*x*z + 4*y*z^2*y - z*x*y + z^2;
>> g = 2*x + 2*x^2 + x*y + 2*y + y*x + y^2;
>> h = x^2 + x*y + y*x + 2*y^2 + z^2 - 1;


% -------------------------------------------------------------------------
% We can check whether the polynomial is a sum of hermitian squares and
% compute the desired decomposition.
% -------------------------------------------------------------------------

>> [IsSohs,X,base,sohs,fsos] = NCsos(f)

***** NCSOStools: module NCSos started *****

Input polynomial has (max) degree 4 and min degree 2.
Detected 5 monomials in 3 variables.
There are 121 monomials in 3 variables of degree at most 4.
There are 117 monomials in 3 variables of degree at most 4 and at least 2.
After Newton Chip Method keeping 3 monomials.
Number of linear constraints: 6.

Starting SDP solver ...
...
Residual norm: 1.1896e-014

Computing SOHS decomposition ... done.
Found SOHS decomposition with 3 factors.

*************** Polynomial is SOHS ***************

IsSohs = 1

X =  1.0000    0.0000   -1.0000
     0.0000    4.0000   -0.0000
    -1.0000   -0.0000    1.0000

base = 'x*y'
       'z*y'
       'z'
 
sohs =    x*y-z
          2*z*y
       5e-008*z
 
fsos = y*x^2*y-y*x*z+4*y*z^2*y-z*x*y+1*z^2

 
% -------------------------------------------------------------------------
% NCsos has many optional parameters. For example we can set the smallest
% value that is considered to be nonzero in numerical calculations.
% -------------------------------------------------------------------------

>> params.precision = 1e-4;
>> [IsSohs,X,base,sohs,fsos] = NCsos(f,params)

***** NCSOStools: module NCSos started *****

Input polynomial has (max) degree 4 and min degree 2.
Detected 5 monomials in 3 variables.
There are 121 monomials in 3 variables of degree at most 4.
There are 117 monomials in 3 variables of degree at most 4 and at least 2.
After Newton Chip Method keeping 3 monomials.
Number of linear constraints: 6.

Starting SDP solver ...
...
Residual norm: 1.1896e-014

Computing SOHS decomposition ... done.
Found SOHS decomposition with 2 factors.

*************** Polynomial is SOHS ***************

IsSohs = 1

X =  1.0000    0.0000   -1.0000
     0.0000    4.0000   -0.0000
    -1.0000   -0.0000    1.0000

base = 'x*y'
       'z*y'
       'z'
 
sohs = x*y-z
       2*z*y
 
fsos = y*x^2*y-y*x*z+4*y*z^2*y-z*x*y+z^2

 
% -------------------------------------------------------------------------
% We can also computes the maximal epsilon such that the polynomial
% g-epsilon is a sum of hermitian squares.
% -------------------------------------------------------------------------

>> g
 
g = 2*x+2*x^2+x*y+2*y+y*x+y^2
 
>> [epsilon,X,base,sohs,gmin] = NCmin(g)

Input polynomial has no constant term ... we add it for computing the base of monomials.
Input polynomial has (max) degree 2 and min degree 0.
Detected 7 monomials in 2 variables.
There are 7 monomials in 2 variables of degree at most 2.
After Newton Chip Method keeping 3 monomials.
Number of linear constraints: 5.

Starting SDP solver ...
...
Residual norm: 2.7182e-010

Computing SOHS decomposition ... done.
Found SOHS decomposition with 3 factors.

epsilon = -1.0000


X = 1.0000    1.0000    1.0000
    1.0000    2.0000    1.0000
    1.0000    1.0000    1.0000

base = ''
       'x'
       'y'

sohs =      1+x+y
                x
       4.1e-007*y
 
gmin = 1+2*x+2*x^2+x*y+2*y+y*x+1*y^2


% -------------------------------------------------------------------------
% NCmin has many optional parameters. For example we can set the smallest
% value that is considered to be nonzero in numerical calculations.
% -------------------------------------------------------------------------

>> params.precision = 1e-4;
>> [epsilon,X,base,sohs,gmin] = NCmin(g,params)

Input polynomial has no constant term ... we add it for computing the base of monomials.
Input polynomial has (max) degree 2 and min degree 0.
Detected 7 monomials in 2 variables.
There are 7 monomials in 2 variables of degree at most 2.
After Newton Chip Method keeping 3 monomials.
Number of linear constraints: 5.

Starting SDP solver ...
...
Residual norm: 2.7182e-010

Computing SOHS decomposition ... done.
Found SOHS decomposition with 2 factors.

epsilon = -1.0000

X = 1.0000    1.0000    1.0000
    1.0000    2.0000    1.0000
    1.0000    1.0000    1.0000

base = ''
       'x'
       'y'
 
sohs = 1+x+y
           x
 
gmin = 1+2*x+2*x^2+x*y+2*y+y*x+y^2


% -------------------------------------------------------------------------
% For given polynomials f and g we can compute the maximal epsilon such
% that the polynomial f-epsilon*g is a sum of hermitian squares.
% -------------------------------------------------------------------------

>> [epsilon,X,base,sohs,poly] = NCdiff(x^2-x^2*y-y*x^2+y*x^2*y+4*y^2,x*y+y*x-2*y*x*y)

Input polynomial has (max) degree 4 and min degree 2.
Detected 8 monomials in 2 variables.
There are 31 monomials in 2 variables of degree at most 4.
There are 28 monomials in 2 variables of degree at most 4 and at least 2.
After Newton Chip Method keeping 3 monomials.
Number of linear constraints: 6.

Starting SDP solver ...
...
Residual norm: 1.285e-011

Computing SOHS decomposition ... done.
Found SOHS decomposition with 3 factors.

***** Optimal solution found:  epsilon_max = 2.000000 *****

epsilon = 2.0000

X =  1.0000   -2.0000   -1.0000
    -2.0000    4.0000    2.0000
    -1.0000    2.0000    1.0000

base = 'x'
       'y'
       'x*y'
 
sohs =                 x-x*y-2*y
       3.95e-006*x*y+5.58e-006*y
                   2.54e-006*x*y
 
poly = x^2-x^2*y-2*x*y-2*y*x-y*x^2+1*y*x^2*y+4*y*x*y+4*y^2


% -------------------------------------------------------------------------
% We can check if the polynomial is convex, i.e., if its second directional
% derivative is a sum of hermitian squares.
% -------------------------------------------------------------------------

>> g
 
g = 2*x+2*x^2+x*y+2*y+y*x+y^2
 
>> [isConvex,g,sohs] = NCisConvex0(g)

isConvex = 1

g = 4*h1^2+2*h1*h2+2*h2*h1+2*h2^2
 
sohs = 2*h1+h2
            h2


% -------------------------------------------------------------------------
% We can also check if the polynomial is convex without SDP solvers
% -------------------------------------------------------------------------

>> isConvex = NCisConvex(g)

isConvex = 1