Decomposition algorithm
using TensorDec
using LinearAlgebra
X = @ring x1 x2 x33-element Array{DynamicPolynomials.PolyVar{true},1}:
x1
x2
x3We 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.0L1 = 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.0The rank of $H_{\sigma}$ will give us an idea on the dimension of $\mathcal{A}_\sigma$.
rank(H)3We 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
x2H0 = 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.0rank(H0)3Let 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.0M = [ 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.14286The 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.635001The 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.0D[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.0D[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.0Looking 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.0We 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'*B03-element Array{DynamicPolynomials.Polynomial{true,Float64},1}:
-0.14285714285714324x1 - 0.2857142857142862x2 + 1.142857142857143
0.28571428571428614x1 - 0.4285714285714279x2 + 0.7142857142857121
-0.14285714285714332x1 + 0.7142857142857134x2 - 0.8571428571428543We deduce the weights $w_i=\sigma(u_i)$:
w = hankel(sigma, U, [L1[1]])3×1 Array{Float64,2}:
1.999999999999992
-1.0000000000000018
5.000000000000002Using 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])Xi3×3 Array{Float64,2}:
4.0 2.0 -1.0
2.0 3.0 1.0
1.0 1.0 1.0w3-element Array{Float64,1}:
-1.0000000000000129
5.000000000000011
1.9999999999999998The series decomposes as $2 \mathfrak{e}_{(-1,1,1)} + 5 \mathfrak{e}_{(2,3,1)} - \mathfrak{e}_{(4,2,1)}$.