Low rank approximation of symmetric tensors

Let us take a random symmetric tensor normally distributed with complex coefficients of order 4 and dimension 3 (the generic symmetric rank is 5), and let us compute by the Riemannian Newton algorithm "rnentr" and the Riemannian Gauss–Newton algorithm "rgnvtr" initialized by a random initial point obeying normal distribution an approximated rank-3 symmetric tensor.

using TensorDec
using DynamicPolynomials
using LinearAlgebra # Take a random symmetric tensor
using Tensors
n = 3; d = 4; r = 3
T = randn(SymmetricTensor{d, n})+randn(SymmetricTensor{d, n})*im
T = convert(Array,T)
# show the first 3 arrays of T
T[:,:,:,1]

3×3×3 Array{ComplexF64, 3}:
[:, :, 1] =
 -0.907877-0.704053im   1.63111+1.9166im    0.707683-1.10415im
   1.63111+1.9166im    0.356316-0.864636im  0.266712-1.17115im
  0.707683-1.10415im   0.266712-1.17115im   -1.00707+1.11311im

[:, :, 2] =
     0.888359+1.18551im   -0.000287941+0.994591im   0.063549+0.920767im
 -0.000287941+0.994591im       0.44311+1.95581im   -0.571079+1.21967im
     0.063549+0.920767im     -0.571079+1.21967im    0.478346-1.00253im

[:, :, 3] =
 0.474927+0.225965im   0.345004+0.096659im  0.378529+0.621547im
 0.345004+0.096659im  -0.454708-0.323135im  0.326064-0.518541im
 0.378529+0.621547im   0.326064-0.518541im  0.988977+1.29563im

# Take the associate homogeneous polynomial P to T by applying the function ahp (for associate homogeneous polynomial)
X = (@polyvar x[1:n])[1]
P = ahp(T, X);

# Take an initial point 
w = ones(r) + fill(0.0+0.0im,r);
V = randn(ComplexF64,n,r);

Applying Riemannian Newton Exact (rne_n_tr) iterations:

w_end, V_end, InfoRNE = approximate(P, w, V; iter = :RNE)

(ComplexF64[5.05399865445625 + 0.0im, 2.3774909931035793 + 0.0im, 5.564701632935593 + 0.0im], ComplexF64[-0.25708654259621544 - 0.20934184778311407im 0.26528717387960404 + 0.16372231300068174im 0.13244665877211262 + 0.2861924486462175im; 0.23097764100261176 - 0.1401358654004084im 0.6261850699359142 - 0.48265493517648883im 0.06264828073916878 + 0.4549788626666184im; 0.655360519402606 + 0.6225723721595502im 0.008275536105741 - 0.5269589232966605im -0.8262225174137973 + 0.0835316950828236im], Dict{String, Real}("d*" => 7.220881621837435, "d0" => 40.4276760623482, "nIter" => 16, "epsIter" => 0.001, "maxIter" => 500))

The reported error is the apolar norm between P and the approximated polynomial.

The initial error is:

InfoRNE["d0"]

40.4276760623482

The final error is:

InfoRNE["d*"]

7.220881621837435

The number of iterations is:

InfoRNE["nIter"]

16

Applying Riemannian Gauss Newton (rgn_v_tr) iterations, one get

w_end, V_end, InfoRGN = approximate(P, w, V; iter = :RGN)

(ComplexF64[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im], ComplexF64[-0.5659979695826164 + 0.10726026844883924im 0.5537406971164952 - 0.34006502788737575im -0.6679966971954102 - 0.05400208838285853im; -0.027432255695650554 - 0.8758998721814316im -0.7730184477273698 - 0.14900174020833434im -0.23146942680034654 + 0.243865447597017im; -0.6287959782796702 + 0.7614820898258986im -0.6721479704250028 - 0.8450340224335604im -0.044370920631584404 - 1.2234092058320365im], Dict{String, Real}("d*" => 4.379474567662476, "d0" => 40.4276760623482, "nIter" => 218, "epsIter" => 0.001, "maxIter" => 500))

The reported error is the apolar norm between P and the approximated polynomial.

The initial error is:

InfoRGN["d0"]

40.4276760623482

The final error is:

InfoRGN["d*"]

4.379474567662476

The number of iterations is:

InfoRGN["nIter"]

218