# improve_factorization.mpl # # By: John May # # Created: March 15th, 2005 # Last Revision: March 18th, 2005 # Method to compute the jacobian of the multiplication map # There are probably faster ways to implement this, but speed doesn't seem to be an issue. multiplication_map_jacobian := proc(fs::list(polynom), # List of Bivariate polynomials vars::list(name) # the variables of fs as a list ) local ns, # degrees of fs numfacs, seqb; # This method uses total degree support of the polynomials # Computes the jacobian of the map (fi)->prod(f_i) evaluated at the given (f_i) numfacs := nops(fs);#print('numfacs'=numfacs); ns := [seq(degree(fs[i],vars),i=1..numfacs)];#print('degrees'=ns); # works for two polys #Matrix([multiplication_map(fs[1], ns[2], vars), multiplication_map(fs[2], ns[1], vars)]); seqb := (n,m) -> [seq(i,i=1..n-1),seq(i,i=n+1..m)]; Matrix([seq(multiplication_map(convert([seq(fs[i],i=seqb(j,numfacs))],`*`), ns[j], vars),j=1..numfacs)]); end proc; # Old multiplication_map_jacobian_old := proc(fs::list(polynom), # List of Bivariate polynomials vars::list(name) # the variables of fs as a list ) local Fs, # Symbolic version of fs ns, # degrees of fs F, # The coeff vector of the symbolic product of the Fs C, # The Jacobian of F eqnvars, numfacs; # This method uses total degree support of the polynomials # Computes the jacobian of the map (fi)->prod(f_i) evaluated at the given (f_i) numfacs := nops(fs);#print('numfacs'=numfacs); ns := [seq(degree(fs[i],vars),i=1..numfacs)];#print('degrees'=ns); # Build symbolic polynomials, so we can take the derivative. # This line could be modified to support > 2 vars Fs := [seq(add(add(a[k,j,i-j]*vars[1]^j*vars[2]^(i-j),j=0..i),i=0..ns[k]), k=1..numfacs)]; # Build the list of variables in the right order eqnvars := convert(Vector([seq(bv2veco(Fs[i],vars),i=1..numfacs)]), list); # Build multiplication map in terms of symbolic coefficients F := Vector([bv2veco(expand(product(Fs[i],i=1..numfacs)),vars)]); C := VectorCalculus[Jacobian](F, eqnvars); # Specialize the Jacobian to the inputs # This line could be modifed to support >2 vars subs({ seq(seq(seq(a[k,j,i-j] = coeff(coeff(fs[k],vars[1],j),vars[2],i-j),j=0..i), i=0..ns[k]),k=1..numfacs)},C); end proc: # multiplication_map_jacobian improve_factorization := proc( f::polynom, # Bivariate polynomial being factored facts::list(polynom), # a list of approximate factors of f vars::list(name) # the variables of f and facs ) local c, n, beo, l, facs, ii, zn, kk, jj,UU,M_M,_SS,_VV,num_S,i,M_new_inverse; # Backward error of the input factorization: c,beo:=backward_error(f, convert(facts,`*`),vars); userinfo(1, {appfac,imfac}, print('org_error'=beo)); userinfo(3, {appfac,imfac}, print('c'=c)); n := nops(facts); # compute the lengths of the coeff. vectors of the factors # this is a stupid way to do this - there is a formula to compute this l:=[seq(op(1,bv2veco(facts[i], vars)),i=1..n)]; userinfo(4, {appfac,imfac}, print('l' = l)); Digits := floor(evalhf(Digits)); # Process the factors: remove extraneous imaginary parts facs:=map(z1->vec2bvo(map(z2->`if`(Im(z2)=0,Re(z2),z2), bv2veco(z1,vars)), vars),facts): # multiply each factor by part of the optimal scaling computed by backward_error() # possibly better ways to do this. if c < 0 then c := -c; facs[1]:=-1.*facs[1] end if; facs:=map(z->evalf(c^(1/n)) * z, facs): # Gauss-Newton Iteration to refine the factors for ii from 1 to 15 do userinfo(4, {appfac,imfac}, print(ii,'bwerror'=norm(expand(f-convert(facs,`*`)),2))) ; # Using using LeastSquares should be more stable, but doesn't seem to work at all. # Zyang modify on Mar.18,2007. # Zyang replaced matrix inverse computation with SVD. Since the matrix is near singular for some special example. # In some case, SVD does not work. See example 8. Hence, I still use MatrixInverse. Apr.18,2007 zn := LinearAlgebra:-MatrixInverse( multiplication_map_jacobian(facs,vars)) . bv2veco(expand(convert(facs,`*`) - f),vars); #M_M:=multiplication_map_jacobian(facs,vars); #UU,_SS,_VV:=LinearAlgebra:-SingularValues(M_M,output=['U','S', 'Vt']); #num_S:=0; #for i from 1 to LinearAlgebra:-ColumnDimension(M_M) do # if _SS[LinearAlgebra:-ColumnDimension(M_M)-i+1]1/x,LinearAlgebra:-SubVector(_SS, # 1..LinearAlgebra:-ColumnDimension(M_M)-num_S)),Vector(num_S)]); #_SS:=LinearAlgebra:-DiagonalMatrix(_SS,LinearAlgebra:-RowDimension(M_M),LinearAlgebra:-RowDimension(M_M)); # M_new_inverse:=LinearAlgebra:-Transpose(_VV).LinearAlgebra:-SubMatrix(_SS.LinearAlgebra:-Transpose(UU),1..LinearAlgebra:-ColumnDimension(M_M),1..-1); #zn:=M_new_inverse.bv2veco(expand(convert(facs,`*`) - f),vars); #zn := LinearAlgebra:-LeastSquares( multiplication_map_jacobian(facs,vars) , bv2veco(expand(convert(facs,`*`) - f),vars)); # Zyang modify on Mar.18, 2007 userinfo(4, {appfac,imfac}, print(ii,'zn'=LinearAlgebra:-Norm(zn,2))); kk:=0; for jj from 1 to n do facs[jj]:=facs[jj]-vec2bvo(zn[ (kk+1)..(kk+l[jj]) ], vars); kk:=kk+l[jj]; end do; if(LinearAlgebra:-Norm(zn,2) < 5*10^(-12) ) then break end if; end do: userinfo(1, {appfac,imfac}, print('iters'=ii)); userinfo(1, {appfac,imfac}, print('new_error'= backward_error(f, convert(facs,`*`), vars)[2], 'improvement' = beo / backward_error(f, convert(facs,`*`), vars)[2])); return facs; end proc: # improve_factorization ################################################################################################################## ## Zyang add procedure "multiplication_map" to this file on Apr.17,2007. ################################################################################################################## bvcauchy := proc(f, deg) local i, j, td, unknowns, g, eqns, __v; # Protect the unknowns td := degree(f, {x,y}); unknowns:=[]: g:=0; for i from deg to 0 by -1 do for j from 0 to i do g := g+__v[j,i-j]*x^j*y^(i-j); unknowns := [op(unknowns), __v[j,i-j]]; end do; end do; eqns := expand(f*g); eqns := convert(bv2vec(eqns),list); eqns := map(x->x=0, eqns); LinearAlgebra[GenerateMatrix](eqns, unknowns)[1]; end proc; multiplication_map := proc(f::polynom, deg::integer, vars::list(name)) local tmpcoeff, td, fterms, gterms, dims, M, rc, cc, d, c, i, j; tmpcoeff:=map(PolynomialTools[CoefficientList], PolynomialTools[CoefficientList](f,vars[1]), vars[2]); td := degree(f, convert(vars,list)); fterms := (td+1)*(td+2)/2; gterms := (deg+1)*(deg+2)/2; dims := [fterms+deg*td+gterms-1, gterms]; M := Matrix(op(dims)); # Build the matrix starting at the last column # Structure: # Each column corresponds to a term of g with constant term last # For a column assoc. with tdeg d term of g the non-zero entries are: # constant term of f + d zeros + tdeg 1 terms of f + d zeros + .. # The bottom lower triangle is zero rc := dims[1]; # row counter cc := dims[2]; # column counter for d from 0 to deg do # do all the columns of total degree d for c from 0 to d do # there are d of these columns # Insert the coeffs of f with the right gaps for i from 0 to td do # terms of tdeg = i for j from i by -1 to 0 do # Insert coeff of x^j*y^(i-j) M[rc, cc] := tmpcoeff[j+1][i-j+1]; rc := rc - 1; end do; # j if i <> td then rc := rc - d; # add a gap end if; end do; # i cc := cc - 1; rc := dims[1] - (dims[2] - cc); end do; # c end do; # d return M; end proc;