Decomposition algorithm
using MultivariateSeries, LinearAlgebra
X = @ring x1 x2 x33-element Vector{DynamicPolynomials.PolyVar{true}}:
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)-26.0dx1^3 + 30.0dx1^2dx2 + 6.0dx1^2dx3 + 72.0dx1*dx2^2 + 20.0dx1*dx2*dx3 + 4.0dx1*dx3^2 + 129.0dx2^3 + 43.0dx2^2dx3 + 15.0dx2*dx3^2 + 6.0dx3^3 + 6.0dx1^2 + 20.0dx1*dx2 + 4.0dx1*dx3 + 43.0dx2^2 + 15.0dx2*dx3 + 6.0dx3^2 + 4.0dx1 + 15.0dx2 + 6.0dx3 + 6.0L1 = monomials(X,0:1)
L2 = monomials(X,0:2)10-element Vector{DynamicPolynomials.Monomial{true}}:
1
x1
x2
x3
x1²
x1x2
x1x3
x2²
x2x3
x3²H = hankel(sigma,L1,L2)4×10 Matrix{Float64}:
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 Vector{DynamicPolynomials.Monomial{true}}:
1
x1
x2H0 = hankel(sigma, B0, B0)3×3 Matrix{Float64}:
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 Matrix{Float64}:
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 Matrix{Float64}:
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, Matrix{Float64}, Vector{Float64}}
values:
3-element Vector{Float64}:
-0.9999999999999991
2.000000000000002
4.000000000000002
vectors:
3×3 Matrix{Float64}:
0.963087 -0.762001 -0.811107
-0.120386 -0.127 -0.324443
-0.240772 0.635001 0.486664The matrices $M_{x_i}$ are diagonal in this basis:
D = [E^(-1)*M[i]*E for i in 1:3]
D[1]3×3 Matrix{Float64}:
-1.0 -3.44169e-15 -7.10543e-15
-4.88498e-15 2.0 -4.44089e-15
4.66294e-15 -4.44089e-15 4.0D[2]3×3 Matrix{Float64}:
1.0 -1.22125e-15 -6.66134e-16
-3.55271e-15 3.0 3.55271e-15
1.11022e-15 2.22045e-15 2.0D[3]3×3 Matrix{Float64}:
1.0 2.22045e-16 2.22045e-16
6.66134e-16 1.0 -8.88178e-16
-9.4369e-16 8.88178e-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 Matrix{Float64}:
-1.0 2.0 4.0
1.0 3.0 2.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 Vector{DynamicPolynomials.Polynomial{true, Float64}}:
-0.1428571428571432x1 - 0.2857142857142861x2 + 1.1428571428571428
-0.14285714285714343x1 + 0.7142857142857132x2 - 0.8571428571428539
0.2857142857142862x1 - 0.4285714285714282x2 + 0.7142857142857126We deduce the weights $w_i=\sigma(u_i)$:
w = hankel(sigma, U, [L1[1]])3×1 Matrix{Float64}:
1.999999999999993
5.000000000000002
-1.0000000000000027Using the command decompose, we can get directly the same decomposition:
w, Xi = decompose(sigma)([1.9999999999999971, 5.0000000000000115, -0.9999999999999997], [-0.9999999999999993 2.0000000000000004 3.999999999999999; 1.0000000000000002 3.000000000000001 1.9999999999999984; 1.0 1.0 1.0])Xi3×3 Matrix{Float64}:
-1.0 2.0 4.0
1.0 3.0 2.0
1.0 1.0 1.0w3-element Vector{Float64}:
1.9999999999999971
5.0000000000000115
-0.9999999999999997The series decomposes as $2 \mathfrak{e}_{(-1,1,1)} + 5 \mathfrak{e}_{(2,3,1)} - \mathfrak{e}_{(4,2,1)}$.