GitHub

Best rank one approximation and optimization on the sphere

using TensorDec
using DynamicPolynomials

# Define the parameters
X = @polyvar x1 x2 x3

# P is a homogeneous polynomial of degree 4 in 3 variables
P = (x1+x2+0.75*x3)^4+1.5*(x1-x2)^4-2*(x1-x3)^4

$

0.5x1^{4} - 2.0x1^{3}x2 + 11.0x1^{3}x3 + 15.0x1^{2}x2^{2} + 9.0x1^{2}x2x3 - 8.625x1^{2}x3^{2} - 2.0x1x2^{3} + 9.0x1x2^{2}x3 + 6.75x1x2x3^{2} + 9.6875x1x3^{3} + 2.5x2^{4} + 3.0x2^{3}x3 + 3.375x2^{2}x3^{2} + 1.6875x2x3^{3} - 1.68359375x3^{4} $

The graph of P in polar coordinates on the sphere looks like this:

title

Let us compute a rank-1 approximation of P. We will compute an initial point by the method SMD, the Julia function that corresponds to this method in the package TensorDec is called "decompose", then we will use the Riemannian Newton algorithm with trust region scheme for the real case, the corresponding Julia function in the package TensorDec is called "rnentr_r".

# Compute an initial point
w1, V1 = decompose(P,1)
([6.545493350772277], [-0.6072039195594636; -0.6330986129014177; -0.4800932684530426;;], Dict{String, Any}("diagonalization" => Dict{String, Any}("case" => "1x1")))

Let us refine this point by using a few number of iterations of rnentr_r, for example 5 iterations.

w_end, V_end, Info = approximate(P, w1, V1; iter= :RNE)
([6.565249939277385], [-0.6236914922916466; -0.624802695988649; -0.4697132247747928;;], Dict{String, Real}("d*" => 9.695336978192456, "d0" => 9.699451819228086, "nIter" => 4, "epsIter" => 0.001, "maxIter" => 500))

The weight in absolute value given by rnentr_r initialized by decompose for rank-1 symmetric tensor approximation is: 7.701576525649196.

The unit vector given by rnentr_r initialized by decompose for rank-1 symmetric tensor approximation is: [0.68061769889553; 0.04831896554028091; -0.7310436550024019].

Comparing with polynomial optimization

Let us use now polynomial optimization to get the spectral norm of P. We use the package MomentToolsand the SDP solver CSDP. To compute this specttral norm, we minimize and maximize P on the unit sphere. The maximum evaluation of P in absolute value on the unit sphere (and that is why we have to use both maximize and minimize functions) gives the spectral norm of P and equivalently a best rank-1 approximation of the symmetric tensor associated to P.

using MomentTools
using CSDP, JuMP
# The function "Optimizer" is a global optimization solver based on positive semi-definite programming
optimizer = CSDP.Optimizer
#using LinearAlgebra

v1, M1 = minimize(P, [x1^2+x2^2+x3^2-1], [], X, 8, optimizer);
v2, M2 = maximize(P, [x1^2+x2^2+x3^2-1], [], X, 8, optimizer);
CSDP 6.2.0
Iter:  0 Ap: 0.00e+00 Pobj:  0.0000000e+00 Ad: 0.00e+00 Dobj:  0.0000000e+00 
Iter:  1 Ap: 6.27e-01 Pobj: -9.1377312e+00 Ad: 6.65e-01 Dobj:  6.5822601e+00 
Iter:  2 Ap: 1.00e+00 Pobj: -8.9378911e+01 Ad: 5.42e-01 Dobj:  6.7430765e+00 
Iter:  3 Ap: 1.00e+00 Pobj: -7.8427861e+01 Ad: 8.83e-01 Dobj:  1.1629773e+00 
Iter:  4 Ap: 1.00e+00 Pobj: -4.7706795e+01 Ad: 8.97e-01 Dobj:  3.8603631e-01 
Iter:  5 Ap: 1.00e+00 Pobj: -1.8751219e+01 Ad: 7.72e-01 Dobj:  5.1044080e-02 
Iter:  6 Ap: 7.86e-01 Pobj: -8.9568745e+00 Ad: 9.24e-01 Dobj: -1.6229630e+00 
Iter:  7 Ap: 9.83e-01 Pobj: -8.6192083e+00 Ad: 7.44e-01 Dobj: -6.4997642e+00 
Iter:  8 Ap: 7.68e-01 Pobj: -7.8344478e+00 Ad: 7.62e-01 Dobj: -7.2868769e+00 
Iter:  9 Ap: 9.65e-01 Pobj: -7.7124187e+00 Ad: 8.41e-01 Dobj: -7.6028422e+00 
Iter: 10 Ap: 9.33e-01 Pobj: -7.7024600e+00 Ad: 9.51e-01 Dobj: -7.6931385e+00 
Iter: 11 Ap: 9.37e-01 Pobj: -7.7016436e+00 Ad: 9.79e-01 Dobj: -7.7012173e+00 
Iter: 12 Ap: 9.65e-01 Pobj: -7.7015833e+00 Ad: 1.00e+00 Dobj: -7.7015626e+00 
Iter: 13 Ap: 9.39e-01 Pobj: -7.7015799e+00 Ad: 1.00e+00 Dobj: -7.7015768e+00 
Iter: 14 Ap: 1.00e+00 Pobj: -7.7015796e+00 Ad: 1.00e+00 Dobj: -7.7015805e+00 
Iter: 15 Ap: 1.00e+00 Pobj: -7.7015795e+00 Ad: 1.00e+00 Dobj: -7.7015804e+00 
Iter: 16 Ap: 1.00e+00 Pobj: -7.7015795e+00 Ad: 1.00e+00 Dobj: -7.7015796e+00 
Iter: 17 Ap: 1.00e+00 Pobj: -7.7015795e+00 Ad: 1.00e+00 Dobj: -7.7015795e+00 
Iter: 18 Ap: 9.60e-01 Pobj: -7.7015795e+00 Ad: 9.55e-01 Dobj: -7.7015795e+00 
Success: SDP solved
Primal objective value: -7.7015795e+00 
Dual objective value: -7.7015795e+00 
Relative primal infeasibility: 4.06e-14 
Relative dual infeasibility: 1.75e-10 
Real Relative Gap: -2.67e-11 
XZ Relative Gap: 1.00e-09 
DIMACS error measures: 7.27e-14 0.00e+00 2.69e-10 0.00e+00 -2.67e-11 1.00e-09
CSDP 6.2.0
Iter:  0 Ap: 0.00e+00 Pobj:  0.0000000e+00 Ad: 0.00e+00 Dobj:  0.0000000e+00 
Iter:  1 Ap: 6.08e-01 Pobj: -8.9542595e+00 Ad: 6.64e-01 Dobj: -2.1488584e+01 
Iter:  2 Ap: 1.00e+00 Pobj: -9.4142852e+01 Ad: 5.33e-01 Dobj: -1.5637445e+01 
Iter:  3 Ap: 1.00e+00 Pobj: -8.1334623e+01 Ad: 8.84e-01 Dobj: -2.2425292e+00 
Iter:  4 Ap: 1.00e+00 Pobj: -5.0436110e+01 Ad: 9.00e-01 Dobj: -1.2447457e+00 
Iter:  5 Ap: 1.00e+00 Pobj: -2.0096267e+01 Ad: 7.90e-01 Dobj: -1.4739780e+00 
Iter:  6 Ap: 8.12e-01 Pobj: -9.3624346e+00 Ad: 8.98e-01 Dobj: -2.4116592e+00 
Iter:  7 Ap: 6.62e-01 Pobj: -6.7261120e+00 Ad: 8.43e-01 Dobj: -4.2835821e+00 
Iter:  8 Ap: 1.00e+00 Pobj: -7.2082331e+00 Ad: 5.32e-01 Dobj: -5.8373968e+00 
Iter:  9 Ap: 9.20e-01 Pobj: -6.6141990e+00 Ad: 9.19e-01 Dobj: -6.4775431e+00 
Iter: 10 Ap: 9.53e-01 Pobj: -6.5676395e+00 Ad: 9.60e-01 Dobj: -6.5604450e+00 
Iter: 11 Ap: 9.57e-01 Pobj: -6.5653593e+00 Ad: 9.82e-01 Dobj: -6.5650429e+00 
Iter: 12 Ap: 9.67e-01 Pobj: -6.5652544e+00 Ad: 9.20e-01 Dobj: -6.5652369e+00 
Iter: 13 Ap: 8.76e-01 Pobj: -6.5652507e+00 Ad: 9.85e-01 Dobj: -6.5652503e+00 
Iter: 14 Ap: 5.02e-01 Pobj: -6.5652503e+00 Ad: 1.00e+00 Dobj: -6.5652504e+00 
Iter: 15 Ap: 1.00e+00 Pobj: -6.5652500e+00 Ad: 1.00e+00 Dobj: -6.5652501e+00 
Iter: 16 Ap: 1.00e+00 Pobj: -6.5652500e+00 Ad: 1.00e+00 Dobj: -6.5652501e+00 
Iter: 17 Ap: 1.00e+00 Pobj: -6.5652500e+00 Ad: 1.00e+00 Dobj: -6.5652500e+00 
Iter: 18 Ap: 9.60e-01 Pobj: -6.5652500e+00 Ad: 9.44e-01 Dobj: -6.5652500e+00 
Success: SDP solved
Primal objective value: -6.5652500e+00 
Dual objective value: -6.5652500e+00 
Relative primal infeasibility: 6.88e-15 
Relative dual infeasibility: 1.23e-09 
Real Relative Gap: -2.17e-10 
XZ Relative Gap: 9.80e-09 
DIMACS error measures: 1.23e-14 0.00e+00 1.88e-09 0.00e+00 -2.17e-10 9.80e-09

The minimum evaluation of P on the unit sphere is: -7.701579459519532.

The maximum evaluation of P on the unit sphere is: 6.565249952183416.

Thus, the maximum weight in absolute value, which is the spectral norm of P, is: 7.701579459519532.

The unit vectors that give this value are: [0.6805571886747267, 0.048429787707634814, -0.7310926539221407] and [-0.6805571886747267, -0.048429787707634814, 0.7310926539221407].

We verify that the second vector and the second weight correspond to the best rank-1 approximation $w_{end}(v_{end}^tX)^4$ of P, given by "rnentr_r" iterations.