We propose an iterative algorithm for the numerical computation of sums of
squares of polynomials approximating given data at prescribed interpolation
points. The method is based on the definition of a convex functional G
arising from the dualization of a quadratic regression over the Cholesky
factors of the sum of squares decomposition. In order to justify the
construction, the domain of G, the boundary of the domain and the behavior at
infinity are analyzed in details. When the data interpolate a positive
univariate polynomial, we show that in the context of the Lukacs sum of squares
representation, G is coercive and strictly convex which yields a unique
critical point and a corresponding decomposition in sum of squares. For
multivariate polynomials which admit a decomposition in sum of squares and up
to a small perturbation of size ε, Gε is always
coercive and so it minimum yields an approximate decomposition in sum of
squares. Various unconstrained descent algorithms are proposed to minimize G.
Numerical examples are provided, for univariate and bivariate polynomials