Re: [Maxima-discuss] complete_square

[Barton W. <wi...@un...>]

There is a share package for Bernstein polynomial stuff. I don't know the theory (and never did), but
I think there are connections between the polynomial SOS form and Bernstein expansions.

--Barton
________________________________
From: Mike Valenzuela <mic...@gm...>
Sent: Tuesday, June 2, 2015 3:08:11 PM
To: Barton Willis
Cc: Luigi Marino; max...@li...
Subject: Re: [Maxima-discuss] complete_square

Thank you again Barton!

Yes I should use Cholesky decomposition when the coefficients of a quadratic equation form a PD matrix. Higher order might be possible. For more complicated expressions I (or someone else) might have to write a semidefinite programming (SDP) method. This would help Maxima rewrite expressions like x^4-x^2-2*x+2 as a square in terms of x^2,x,1. Does anyone know of a lisp, Maxima, or otherwise-compatible-with-Maxima SDP program?

However, the SDP method appears to find a different set of solutions efficiently. The spectral method (code at the bottom) will complete the square in linear combinations of the variables. Consider x^2+2*x*y+y^2+2*x+2*y. Using the code at the bottom of this email it produces the following sum of squares:
 (1-sqrt(3))*(y+x-(sqrt(3)+1))^2     (sqrt(3)+1)*(y+x+sqrt(3)-1)^2
--------------------------------- + -------------------------------
         (2*(sqrt(3)+3))                        (2*(3-sqrt(3)))

The spectral method is probably worthless as far as Single Use Expressions (SUE) goes, but still useful in some other sense. It is a sum of squares in a linear combination of the variables. It looks like the SDP would need to know this form in advance; perhaps that is a good though. Despite the algebraic eivects method being slow and how large its results can be, the spectral factoring might be of interest to someone. Also, we could introduce a special case where all the coefficients are numbers and use a numerical eigenvector package in that case.


---------
/* This code is released under the GNU General Public License 3 */
/* http://www.gnu.org/licenses/gpl-3.0.en.html */

load(eigen)$

normalized_vec(v):= v / sqrt(lsum(x^2,x,v))$

Spectral_Factor(sym_matrix):= block([Q, Lambda, n, vects, vals],
/* declare local vars */
local(Q, Lambda, n, vects, vals),

/* Error Checking Missing */

/* The bread and butter for many linear algebra techniques */
[vals, vects]: ueivects(sym_matrix),

/* Allocate memory */
Q: genmatrix(lambda([i,j],0),length(sym_matrix),length(sym_matrix)),
Lambda: genmatrix(lambda([i,j],0),length(sym_matrix),length(sym_matrix)),
n:0,

/* Enter all eivects and eivals */
for i:1 thru length(vects) do
  for j:1 thru length(vects[i]) do block(
    n:n+1,
    Q[n]: vects[i][j],
    Lambda[n][n]: vals[1][i]
  ),

/* Orthogonalize Q in case we had repeated values */
if( length(Q)>1 ) then Q: apply(matrix, map(normalized_vec, gramschmidt(Q))),
Q: transpose(radcan(Q)),
[Q, Lambda])$



/* Now test it */
load(scifac)$
load(sqdnst)$
/* I manually encoded the polynomial x^2+2*x*y+y^2+2*x+2*y */
D: matrix([1,1,1],[1,1,1],[1,1,0])$
[Q, Lambda]:Spectral_Factor(D);
D= radcan( sqrtdenest( Q . Lambda . transpose(Q)));
tmp: matrix([x,y,1]) . Q . sqrt(Lambda);

tmp . transpose(tmp);
/* simplification is difficult to do without destroying the factored form */
factored_expression: (scanmap(combine,sqrtdenest( tmp . transpose(tmp) )));
factored_expression: map(lambda([x],x/2), gcfac(2*factored_expression));



On Tue, Jun 2, 2015 at 4:29 AM, Barton Willis <wi...@un...<mailto:wi...@un...>> wrote:

See also (for example): http://homepages.math.uic.edu/~jan/mcs563/sos.pdf


--Barton