Decomposition algorithm

Decomposition algorithm

using TensorDec
using LinearAlgebra
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)
4.0dx1*dx3 + 15.0dx2 + 6.0dx3 + 20.0dx1*dx2 + 6.0dx3^3 + 43.0dx2^2dx3 - 26.0dx1^3 + 129.0dx2^3 + 30.0dx1^2dx2 + 15.0dx2*dx3 + 20.0dx1*dx2*dx3 + 6.0dx1^2 + 6.0dx3^2 + 4.0dx1 + 43.0dx2^2 + 6.0dx1^2dx3 + 4.0dx1*dx3^2 + 72.0dx1*dx2^2 + 15.0dx2*dx3^2 + 6.0
L1 = monoms(X,1)
L2 = monoms(X,2)
10-element Array{DynamicPolynomials.Monomial{true},1}:
 1   
 x1  
 x2  
 x3  
 x1² 
 x1x2
 x1x3
 x2² 
 x2x3
 x3²
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
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

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

v, E = eigen(M[1])
Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}
eigenvalues:
3-element Array{Float64,1}:
 -0.9999999999999991
  4.000000000000002 
  2.000000000000002 
eigenvectors:
3×3 Array{Float64,2}:
  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.99441e-15  -3.66374e-15
  4.21885e-15   4.0          -4.44089e-15
 -4.66294e-15  -3.9968e-15    2.0
D[2]
3×3 Array{Float64,2}:
  1.0          -4.44089e-16  -1.44329e-15
  8.88178e-16   2.0           2.66454e-15
 -3.55271e-15   2.66454e-15   3.0
D[3]
3×3 Array{Float64,2}:
  1.0           3.33067e-16  1.11022e-16
 -9.4369e-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.14285714285714324x1 - 0.2857142857142862x2 + 1.142857142857143 
 0.28571428571428614x1 - 0.4285714285714279x2 + 0.7142857142857121 
 -0.14285714285714332x1 + 0.7142857142857134x2 - 0.8571428571428543

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

w = hankel(sigma, U, [L1[1]])
3×1 Array{Float64,2}:
  1.999999999999992 
 -1.0000000000000018
  5.000000000000002

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.0000000000000129
  5.000000000000011 
  1.9999999999999998

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