GitHub

Decomposition algorithm

using DynamicPolynomials, MultivariateSeries, LinearAlgebra
X = @polyvar x1 x2 x3
(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 + 6.0dx3 + 15.0dx2 + 4.0dx1 + 6.0dx3^2 + 15.0dx2*dx3 + 43.0dx2^2 + 4.0dx1*dx3 + 20.0dx1*dx2 + 6.0dx1^2 + 6.0dx3^3 + 15.0dx2*dx3^2 + 43.0dx2^2dx3 + 129.0dx2^3 + 4.0dx1*dx3^2 + 20.0dx1*dx2*dx3 + 72.0dx1*dx2^2 + 6.0dx1^2dx3 + 30.0dx1^2dx2 - 26.0dx1^3
L1 = monomials(X,0:1)
L2 = monomials(X,0:2)
10-element MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}}:
 1
 x3
 x2
 x1
 x3²
 x2x3
 x2²
 x1x3
 x1x2
 x1²
H = hankel(sigma,L1,L2)
4×10 Matrix{Float64}:
  6.0   6.0  15.0   4.0   6.0  15.0   43.0   4.0  20.0    6.0
  6.0   6.0  15.0   4.0   6.0  15.0   43.0   4.0  20.0    6.0
 15.0  15.0  43.0  20.0  15.0  43.0  129.0  20.0  72.0   30.0
  4.0   4.0  20.0   6.0   4.0  20.0   72.0   6.0  30.0  -26.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 MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}}:
 1
 x3
 x2
H0 = hankel(sigma, B0, B0)
3×3 Matrix{Float64}:
  6.0   6.0  15.0
  6.0   6.0  15.0
 15.0  15.0  43.0
rank(H0)
2

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 Matrix{Float64}:
  4.0   4.0  20.0
  4.0   4.0  20.0
 20.0  20.0  72.0
M = [ H0^(-1)*H[i] for i in 1:3 ]
M[1]
SingularException(2)



Stacktrace:

  [1] checknonsingular

    @ /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/factorization.jl:19 [inlined]

  [2] checknonsingular

    @ /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/factorization.jl:22 [inlined]

  [3] #lu!#170

    @ /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/lu.jl:82 [inlined]

  [4] lu!

    @ /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/lu.jl:80 [inlined]

  [5] #lu#176

    @ /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/lu.jl:299 [inlined]

  [6] lu (repeats 2 times)

    @ /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/lu.jl:298 [inlined]

  [7] inv(A::Matrix{Float64})

    @ LinearAlgebra /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/dense.jl:917

  [8] literal_pow(#unused#::typeof(^), A::Matrix{Float64}, #unused#::Val{-1})

    @ LinearAlgebra /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/generic.jl:1068

  [9] (::var"#3#4")(i::Int64)

    @ Main ./none:0

 [10] iterate

    @ ./generator.jl:47 [inlined]

 [11] collect(itr::Base.Generator{UnitRange{Int64}, var"#3#4"})

    @ Base ./array.jl:782

 [12] top-level scope

    @ In[10]:1

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])
UndefVarError: `M` not defined



Stacktrace:

 [1] top-level scope

   @ In[11]:1

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

D = [E^(-1)*M[i]*E for i in 1:3]
D[1]
UndefVarError: `E` not defined



Stacktrace:

 [1] (::var"#5#6")(i::Int64)

   @ Main ./none:0

 [2] iterate

   @ ./generator.jl:47 [inlined]

 [3] collect(itr::Base.Generator{UnitRange{Int64}, var"#5#6"})

   @ Base ./array.jl:782

 [4] top-level scope

   @ In[12]:1
D[2]
UndefVarError: `D` not defined



Stacktrace:

 [1] top-level scope

   @ In[13]:1
D[3]
UndefVarError: `D` not defined



Stacktrace:

 [1] top-level scope

   @ In[14]:1

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]
UndefVarError: `D` not defined



Stacktrace:

 [1] (::var"#7#8")(::Tuple{Int64, Int64})

   @ Main ./none:0

 [2] iterate

   @ ./generator.jl:47 [inlined]

 [3] collect(itr::Base.Generator{Base.Iterators.ProductIterator{Tuple{UnitRange{Int64}, UnitRange{Int64}}}, var"#7#8"})

   @ Base ./array.jl:782

 [4] top-level scope

   @ In[15]:1

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
UndefVarError: `E` not defined



Stacktrace:

 [1] top-level scope

   @ In[16]:1

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

w = hankel(sigma, U, [L1[1]])
UndefVarError: `U` not defined



Stacktrace:

 [1] top-level scope

   @ In[17]:1

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

w, Xi = decompose(sigma)
([2.0, 5.000000000000065, -1.0000000000000482], [-1.0 1.9999999999999987 4.000000000000004; 1.0 2.9999999999999982 2.0000000000000013; 1.0 1.0 1.0])
Xi
3×3 Matrix{Float64}:
 -1.0  2.0  4.0
  1.0  3.0  2.0
  1.0  1.0  1.0
w
3-element Vector{Float64}:
  2.0
  5.000000000000065
 -1.0000000000000482

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