Getting the Most from Algebraic Solvers in Mathematica

Adam Strzebonski

Wolfram Research Inc.

Abstract

This talk will give a survey of Mathematica functions related to solving algebraic equations and inequalities. It will also discuss the choice of the most appropriate solvers for various types of problems and the ways of formulating the problems so that they can be solved most efficiently.

Solving algebraic equations and inequalities

• Complete description of the solution set
• Generic description of the solution set
• Solution instances
• The number of solutions
• Quantifier elimination & existence of solutions
• Global optimization
• Local optimization

Complete description of the solution set

• Reduce
• Solve with MaxExtraConditions->All
• GroebnerBasis
• CylindricalDecomposition
• RootIntervals
• NRoots

The solution set represented as a Boolean formula.

Reduce[z^2-x^2y==0,{x,y,z},Reals]

The solution set represented as a list of replacement rules.

Solve[z^2-x^2y==0,z,Reals,MaxExtraConditions->All]

Picture

Show[{Plot3D[{-(Sqrt[(x^2) y]),(Sqrt[(x^2) y])},{x,-1,1},{y,0,1}],Graphics3D[{Blue,Thick,Line[{{0,-1,0},{0,0,0}}]}]},PlotRange->All,ViewPoint->{3,-2,2}]

Intervals isolating roots of univariate polynomials.

Real roots:

RootIntervals[x^9-7x^2+3x+1]

All roots:

RootIntervals[x^9-7x^2+3x+1,Complexes]

Picture

poly=x^9-7x^2+3x+1;
ii=RootIntervals[poly,Complexes][[1]];
f=poly/.x->a+b I;
{ref,imf}=ComplexExpand[{Re[f],Im[f]}];
pl=ContourPlot[{ref,imf},{a,-3.2,3.2},{b,-3.2,3.2}];
SeedRandom[21];
ci={Hue[Random[]],Opacity[0.3],Rectangle[{Re[#[[1]]],Im[#[[1]]]},{Re[#[[2]]],Im[#[[2]]]}]}&/@Drop[ii,3];
ri={Hue[Random[]],Thick,Line[{#,0}&/@#]}&/@Take[ii,3];
rts={Red,PointSize[Large],Point[{Re[#],Im[#]}]}&/@(x/.NSolve[poly==0,x]);
Show[{pl,Graphics[{ci,ri,rts}]}]

ii

Generic description of the solution set

• Solve
• NSolve
• GroebnerBasis with CoefficientDomain->RationalFunctions
• GenericCylindricalDecomposition
• LinearSolve
• Roots

Solve[z^2-x^2y==0,z,Reals]

Picture

The dashed lines indicate the missing parts of the solution set.

Show[{Plot3D[{-(Sqrt[(x^2) y]),(Sqrt[(x^2) y])},{x,-1,1},{y,0,1}],Graphics3D[{Red,Dashed,Thick,Line[{{0,-1,0},{0,1,0}}],Line[{{-1,0,0},{1,0,0}}]}]},PlotRange->All,ViewPoint->{3,-2,2}]

Solution instances

• FindInstance
• FindRoot
• LinearSolve
• SemialgebraicComponentInstances

Find a solution instance.

f=(x^8)-7 (x^6) y+4 (x^6) (y^2)+35 (x^4) (y^3)+6 (x^4) (y^4)-21 (x^2) (y^5)+4 (x^2) (y^6)+(y^7)+(y^8);
FindInstance[f<0,{x,y},Reals]

Find at least one point in each conneted component of the solution set (works for systems of real polynomial equations and inequalities).

SemialgebraicComponentInstances[f<0,{x,y}]

Picture

pl=RegionPlot[f<0,{x,-1,1},{y,-1,1},PlotPoints->100];
p0=Graphics[{Red,PointSize[Large],Point[{x,y}]/.FindInstance[f<0,{x,y},Reals]}];
p1=Graphics[{Green,PointSize[Large],Point[{x,y}]/.SemialgebraicComponentInstances[f<0,{x,y}]}];
Show[{pl,p0,p1}]

The number of solutions

• CountRoots

The number of real roots of a polynomial.

f=Sum[RandomInteger[{-(10^6),(10^6)}] (x^(RandomInteger[{1,(10^6)}])),{10}]+RandomInteger[{-(10^6),(10^6)}]
CountRoots[f,x]

The number of real roots of an exp-log function.

CountRoots[(E^x)-2 x Log[x]-x-(7/4),x]

The number of a holomorphic elementary function in a closed rectangle.

CountRoots[Sin[(E^x)]-x+1,{x,-3-3 I,3+3 I}]

Quantifier elimination & existence of solutions

• Resolve
• Reduce
• Solve
• Eliminate
• SolveAlways
• CylindricalDecomposition
• GroebnerBasis

Eliminate the quantifier.

Resolve[(Subscript[∃, x])((a-1) (x^2)+b x+a+b==0&&(a-b) (x^2)+b x+a<0),Reals]

Eliminate the quantifier and solve the resulting system.

Reduce[(Subscript[∃, x])((a-1) (x^2)+b x+a+b==0&&(a-b) (x^2)+b x+a<0),{a,b},Reals]

Global optimization

• Minimize/Maximize
• MinValue/MaxValue
• ArgMin/ArgMax
• NMinimize/NMaximize
• NMinValue/NMaxValue
• NArgMin/NArgMax
• LinearProgramming

For linear problems computing the minimum value has the same cost as computing both the minimum value and the minimum argument. For non-linear problems there may be a difference. In particular, finding the minimum value subject to strict inequality constraints can be significantly faster than finding both the minimum value and the minimum argument.

MinValue[x+y+z,(x^4)+(y^4)+(z^4)<1&&(x^3)>y z+(1/2),{x,y,z}]//AbsoluteTiming

TimeConstrained[Minimize[x+y+z,(x^4)+(y^4)+(z^4)<1&&(x^3)>y z+(1/2),{x,y,z}],120]

Local optimization

• FindMinimum/FindMaximum
• FindMinValue/FindMaxValue
• FindArgMin/FindArgMax

A local minimum.

f=((x^2)-7) ((x^2)-2) ((x^2)-1);
FindMinimum[f,x]

An approximation of a global minimum.

NMinimize[f,x]

An exact value of a global minimum.

Minimize[f,x]

N[%]

Picture

pl=Plot[f,{x,-3,3}];
lm=FindMinimum[f,x];
gm=NMinimize[f,x];
p0=Graphics[{Red,PointSize[Large],Point[{lm[[2,1,2]],lm[[1]]}]}];
p1=Graphics[{Green,PointSize[Large],Point[{gm[[2,1,2]],gm[[1]]}]}];
Show[{pl,p0,p1}]

Generic vs. complete solutions

Reduce returns complete solutions.

Reduce[a x^2+b x+c==0,x]

Solve gives “generic” solutions, that is solutions that do not require equational conditions on continuous parameters.

Solve[a x^2+b x+c==0,x]

Generic solutions may require that continuous parameters satisfy inequalities.

Solve[a x^2+b x+c==0,x,Reals]

Generic vs. complete solutions II

This system has no generic solutions.

Solve[a x^2==b&&c x^3==d,x]

With MaxExtraConditions→Automatic, Solve returns solutions that require the minimal number of equational conditions on continuous parameters.

Solve[a x^2==b&&c x^3==d,x,MaxExtraConditions->Automatic]

With MaxExtraConditions→All, Solve returns all solutions.

Solve[a x^2==b&&c x^3==d,x,MaxExtraConditions->All]

Full-dimensional components of the solution set cannot be expressed through replacement rules.

Solve[x^2<a,x,Reals]

Reduce can describe arbitrary solution sets using Boolean formulas.

Reduce[x^2<a,x,Reals]

Assumptions on parameters in Solve

In Mathematica 7 and earlier Solve treated equations not containing the variables as “assumptions”.

Solve[a^2==1&&x==a,x] (* Evaluated with Mathematica 7 *)

Problems:
• Different solution sets for logically equivalent systems.
• Not clear how to identify assumptions for arbitrary Boolean combinations of equations and inequalities.
• The notion of “genericity wrt. assumptions” is not well defined.

In Mathematica 8 Solve no longer has the “syntactic assumptions” feature.

Solve[a^2==1&&x==a,x]

With MaxExtraConditions→Automatic, Solve returns the solutions given by Mathematica 7, with the necessary assumptions added to each solution.

Solve[a^2==1&&x==a,x,MaxExtraConditions->Automatic]

ConditionalExpression can be removed using Normal.

Normal[%]

Transcendental equations with symbolic parameters

Solve uses non-equivalent transformations in an attempt to find some solutions fast.

Solve[x^a==2,x]

Reduce finds a complete solution with all necessary conditions on parameters (this may be slow).

(red=Reduce[x^a==2,x])//AbsoluteTiming

Cases[red,x==_,Infinity]//Union

c0=red/.{C[1]->0,Re[a]->re,Im[a]->im,x≠0->True,x==2^(1/a)->True}

RegionPlot[c0,{re,-0.5,0.5},{im,-0.5,0.5}]

f=x^a/.x->2^(1/a)

Plot[Evaluate[f/.a->I t],{t,-0.5,0.5}]

Exact vs. approximate solutions (linear equations)

Integer Coefficients

Approximate solutions of linear equations can often be computed faster than exact solutions.

A=RandomInteger[{-1000,1000},{200,200}];
b=RandomInteger[{-1000,1000},{200}];
vars=x/@Range[200];
(sol1=Solve[A.vars==b,vars]);//AbsoluteTiming

(sol2=NSolve[A.vars==b,vars]);//AbsoluteTiming

Max[Abs[(vars/.sol1)-(vars/.sol2)]]

Non-rational numeric Coefficients

The speed difference is especially large for systems with coefficients containing a large number of distinct non-rational numeric expressions. Here is an extreme case where all coefficients are algebraically independent over ℚ (by Lindemann–Weierstrass theorem).

eqns[n_]:=
With[{a=2^(1/(n(n+1)))},
And@@Table[Sum[Exp[a^((j-1)(n+1)+i)]x[i],{i,n}]==Exp[a^(j(n+1))],{j,n}]]
vars[n_]:=x/@Range[n];
eqns[4]

sol1=Table[AbsoluteTiming[Solve[eqns[n],vars[n]]],{n,10}];
sol2=Table[AbsoluteTiming[NSolve[eqns[n],vars[n]]],{n,10}];

Timings Solve vs. NSolve (logarithmic scale)

ListLogPlot[{First/@sol1,First/@sol2},Joined->True]

LeafCount of exact solutions (logarithmic scale)

ListLogPlot[LeafCount[#[[2]]]&/@sol1,Joined->True]

Correctness of inexact solutions

An approximate method may give a wrong result without a warning...

sol2[[10]]

Even with non-machine working precision error estimates given may be too low.

NSolve[eqns[10],vars[10],WorkingPrecision->20]//AbsoluteTiming

eqns[10]/.%[[2]]/.{x[8]->2,x[10]->3}

With a higher working precision we get a correct result, still > 10000 times faster than the exact computation.

NSolve[eqns[10],vars[10],WorkingPrecision->100]//AbsoluteTiming

eqns[10]/.%[[2]]

Exact vs. approximate solutions (polynomial roots)

Another case in which solutions can be obtained by purely approximate computation is univariate polynomial equations.

Computing solutions using machine precision computations we may get very imprecise results without a warning.

f=Expand[Product[(x-i),{i,30}]];
NSolve[f==1,x]

With non-machine precision the results are given with correct error estimates.

NSolve[N[f,16]==1,x]

If the input is exact, the precision of computations is chosen automatically to get sufficiently precise results.

NSolve[f==1,x,WorkingPrecision->16]

Approximating solutions computed with an exact method yields solutions computed to the requested precision.

N[Solve[f==1,x],16]

Exact vs. approximate input

In Mathematica 8 purely numeric solution methods are implemented only for systems of linear equations and for univariate polynomial equations. For all other systems, or when symbolic parameters are present, (partially) symbolic methods are used.

When no numeric solution methods exist inexact coefficients are rationalized and an exact version of the input is solved.

NSolve[0.25Log[x]+0.75Log[x^2]==1,x]

Rationalizing coefficients in an inexact input may not yield a solvable exact problem ...

NSolve[0.25Log[2,x]+0.75Log[4,x^2]==1,x]

... while the original problem given with exact coefficients is solvable.

NSolve[1/4Log[2,x]+3/4Log[4,x^2]==1,x]

Solvers generally work better with rational rather than floating-point number coeffcients, even when computing approximate solutions. Approximate input to a symbolic method can make the computation harder. The precision of computations can be set using the WorkingPrecision option (present in all solvers). Finite working precision will be used in numeric algorithms, while symbolic algorithms will still be able to use the exact version of the input.

When to use approximate coefficients

• The problem can be solved by purely numeric methods (linear or univariate polynomial equations).
• The coefficients contain a large number of distinct non-rational numeric expressions.
• Only approximate values of coefficients are known, the coefficients are given as non-machine precision numbers with meaningful values of error, and error estimates of solutions are needed.

When to use exact coefficients

• Solving the problem requires symbolic methods.
• The coefficients are rational numbers.

Integer solution instances

General methods for finding integer solutions exist only for systems of linear equations and inequalities, univariate polynomial equations, binary quadratic equations and a few other special classes of equations. When none of the general methods apply Mathematica uses one of the following heuristics.

Modular sieve

This method applies to arbitrary systems of polynomial equations and inequalities with integer coefficients. It works by finding solutions of equations are modulo (the product of) a few small primes. The method gives up if it does not find a solution instance or prove that there are no solutions after trying a certain number of solution candidates. This number can be set using the system option “SieveMaxPoints”. The default value is 10000.

FindInstance[2 x^2+3 y^3+4 z^4==4411,{x,y,z},Integers]

SetSystemOptions["ReduceOptions"->{"SieveMaxPoints"->1000000}];
FindInstance[2 x^2+3 y^3+4 z^4==4411,{x,y,z},Integers]

Exhaustive search

This method applies to arbitrary systems of equations and inequalities with integer coefficients. It works by trying all values of variables within finite ranges. The use of the method depends of the value {m, M} of the “ExhaustiveSearchMaxPoints” system option.

"ExhaustiveSearchMaxPoints"/.("ReduceOptions"/.SystemOptions[])

If the system contains explicit finite range specifications for all variables and the total number of points within the range does not exceed m, exhaustive search is used without trying other methods. Otherwise if the system contains explicit finite range specifications for all variables and the total number of points within the range does not exceed M, exhaustive search is used if all other instance finding methods fail (but before calling Reduce). If the system does not contain explicit finite range specifications for all variables but contains polynomial inequalities that imply finite bounds on some of the variables, FindInstance is called recursively for all values of the variable with the smallest range, provided the number of points within this range does not exceed M.

This runs an exhaustive search for 1≤x≤100 and 1≤y≤100.

FindInstance[Abs[Sin[x]-Cos[y]]<(1/1000)&&1≤x≤100&&1≤y≤100,{x,y},Integers]//AbsoluteTiming

This finds the bounds on x and y implied by the inequalities and calls FindInstance recursively for 1≤x≤100. This is slower than the previous example, even though the number of possible values of {x,y} is smaller.

FindInstance[Abs[Sin[x]-Cos[y]]<(1/1000)&&x≥1&&y≥1&&(x^2)+(y^2)≤10000,{x,y},Integers]//AbsoluteTiming

This calls FindInstance recursively for 1≤x≤100. This is even slower, since each recursive call contains no bounds on y.

FindInstance[Abs[Sin[x]-Cos[y]]<(1/1000)&& 1≤x≤100,{x,y},Integers]//AbsoluteTiming

There are no solutions of this inequality for 1≤x≤100 and 1≤y≤100.

FindInstance[Abs[Sin[x]-Cos[y]]<(1/100000)&&1≤x≤100&&1≤y≤100,{x,y},Integers]//AbsoluteTiming

If we increase the range of x and y, the total number of possible values of {x,y} exceeds M=10000 and exhaustive search is not used. FindInstance instead tries to get a full solution of the system using Reduce, which in this case is very slow.

TimeConstrained[FindInstance[Abs[Sin[x]-Cos[y]]<(1/100000)&&1≤x≤1000&&1≤y≤1000,{x,y},Integers],120]

With M=(10^6) an FindInstance uses exhaustive search and finds a solution instance quite fast.

SetSystemOptions["ReduceOptions"->{"ExhaustiveSearchMaxPoints"->{10^3,10^6}}];

FindInstance[Abs[Sin[x]-Cos[y]]<(1/100000)&&1≤x≤1000&&1≤y≤1000,{x,y},Integers]//AbsoluteTiming