Re: [Maxima-discuss] complete_square

[Mike V. <mic...@gm...>]

Thanks for the tip!

Also after briefly reading about Bernstein polynomials they appear as if
they can be used as a certificate of positivity over the interval [0,1].
Each Bernstein polynomial is non-negative in its variables over the
interval [0,1]. So representing a polynomial as a sum of Bernstein
polynomials with non-negative coefficients implies the original polynomial
is non-negative in the interval [0,1]. I'm not sure how to use Bernstein
polynomials to convert into SoS form, but apparently methods to obtain
certificates of positivity is very active research.



On Tue, Jun 2, 2015 at 2:39 PM, Barton Willis <wi...@un...> wrote:

>  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...> wrote:
>
>>  See also (for example):
>> http://homepages.math.uic.edu/~jan/mcs563/sos.pdf
>>
>>
>>   --Barton
>>
>>
>