Decomposition algorithm

Decomposition algorithm

using TensorDec
X = @ring x1 x2 x3
3-element Array{DynamicPolynomials.PolyVar{true},1}:
 x1
 x2
 x3

We want to find a sparse representation of the following series known up to degree 3:

sigma = dual(6.0 + 4.0*x1 + 15.0*x2 + 6.0*x3 + 6.0*x1^2 + 20.0*x1*x2 + 4.0*x1*x3 + 43.0*x2^2 + 15.0*x2*x3 + 6.0*x3^2 - 26.0*x1^3 + 30.0*x1^2*x2 + 6.0*x1^2*x3 + 72.0*x1*x2^2 + 20.0*x1*x2*x3 + 4.0*x1*x3^2 + 129.0*x2^3 + 43.0*x2^2*x3 + 15.0*x2*x3^2 + 6.0*x3^3)
6.0 + 20.0dx1*dx2*dx3 + 6.0dx1^2 + 15.0dx2*dx3 + 129.0dx2^3 + 43.0dx2^2dx3 + 4.0dx1*dx3 + 6.0dx3^3 + 4.0dx1*dx3^2 + 20.0dx1*dx2 + 30.0dx1^2dx2 + 6.0dx3 + 6.0dx1^2dx3 - 26.0dx1^3 + 15.0dx2*dx3^2 + 4.0dx1 + 15.0dx2 + 43.0dx2^2 + 72.0dx1*dx2^2 + 6.0dx3^2
L1 = monoms(X,1)
L2 = monoms(X,2)
10-element Array{DynamicPolynomials.Monomial{true},1}:
 1   
 x1  
 x2  
 x3  
 x1^2
 x1x2
 x1x3
 x2^2
 x2x3
 x3^2
H = hankel(sigma,L1,L2)
4×10 Array{Float64,2}:
  6.0   4.0  15.0   6.0    6.0  20.0   4.0   43.0  15.0   6.0
  4.0   6.0  20.0   4.0  -26.0  30.0   6.0   72.0  20.0   4.0
 15.0  20.0  43.0  15.0   30.0  72.0  20.0  129.0  43.0  15.0
  6.0   4.0  15.0   6.0    6.0  20.0   4.0   43.0  15.0   6.0

The rank of $H_{\sigma}$ will give us an idea on the dimension of $\mathcal{A}_\sigma$.

rank(H)
3

We check that $\{1, x_1, x_2\}$ is a basis of $\mathcal{A}_\sigma$:

B0 = L1[1:3]
3-element Array{DynamicPolynomials.Monomial{true},1}:
 1 
 x1
 x2
H0 = hankel(sigma, B0, B0)
3×3 Array{Float64,2}:
  6.0   4.0  15.0
  4.0   6.0  20.0
 15.0  20.0  43.0
rank(H0)
3

Let us compute the shifted (truncated) Hankel operators.

H1 = hankel(sigma, B0, B0*x1)
H2 = hankel(sigma, B0, B0*x2)
H3 = hankel(sigma, B0, B0*x3);
H  = [H1,H2,H3]
H[1]
3×3 Array{Float64,2}:
  4.0    6.0  20.0
  6.0  -26.0  30.0
 20.0   30.0  72.0

We deduce the operators of multiplication by the variables in the basis $B_0$:

M = [ H0^(-1)*H[i] for i in 1:3 ]
M[1]
3×3 Array{Float64,2}:
  1.11022e-16   9.14286  -0.571429
  1.0           3.85714   1.57143 
 -1.11022e-16  -4.28571   1.14286

The eigenvalues and eigenvectors of $M_{x_1}$ are

v, E = eig(M[1])
([-1.0, 4.0, 2.0], [0.963087 -0.811107 -0.762001; -0.120386 -0.324443 -0.127; -0.240772 0.486664 0.635001])

The matrices $M_{x_i}$ are diagonal in this basis:

D = [E^(-1)*M[i]*E for i in 1:3]
D[1]
3×3 Array{Float64,2}:
 -1.0          -6.21725e-15  -2.55351e-15
  4.44089e-15   4.0          -4.44089e-15
 -4.21885e-15  -3.9968e-15    2.0
D[2]
3×3 Array{Float64,2}:
  1.0          -4.44089e-16  -6.66134e-16
  9.99201e-16   2.0           2.66454e-15
 -3.10862e-15   2.66454e-15   3.0
D[3]
3×3 Array{Float64,2}:
  1.0           4.44089e-16  4.44089e-16
 -8.32667e-16   1.0          6.66134e-16
  5.55112e-16  -6.66134e-16  1.0

Looking at the corresponding terms on the diagonal, we get the coordinates of the points $\Xi$:

Xi = [ D[i][j,j] for i in 1:3, j in 1:3]
3×3 Array{Float64,2}:
 -1.0  4.0  2.0
  1.0  2.0  3.0
  1.0  1.0  1.0

We normalize the eigenvectors by $v_i \over v_i(\xi_i)$ and get the interpolation polynomials at the points $\xi_i$:

Dg = E'*vcat(fill(1.,1,3), Xi[1:2,:])
E = E*Dg^(-1)
U = E'*B0
3-element Array{DynamicPolynomials.Polynomial{true,Float64},1}:
 -0.142857x1 - 0.285714x2 + 1.14286 
 0.285714x1 - 0.428571x2 + 0.714286 
 -0.142857x1 + 0.714286x2 - 0.857143

We deduce the weights $w_i=\sigma(u_i)$:

w = hankel(sigma, U, [L1[1]])
3×1 Array{Float64,2}:
  2.0
 -1.0
  5.0

Using the command decompose, we can get directly the same decomposition:

w, Xi = decompose(sigma)
([-1.0, 5.0, 2.0], [4.0 2.0 -1.0; 2.0 3.0 1.0; 1.0 1.0 1.0])
Xi
3×3 Array{Float64,2}:
 4.0  2.0  -1.0
 2.0  3.0   1.0
 1.0  1.0   1.0
w
3-element Array{Float64,1}:
 -1.0
  5.0
  2.0

The series decomposes as $2 \mathfrak{e}_{(-1,1,1)} + 5 \mathfrak{e}_{(2,3,1)} - \mathfrak{e}_{(4,2,1)}$.