Find the real roots of polynomials in SAS

0

A SAS programmer had many polynomials for which he wanted to compute the real roots. By the Fundamental Theorem of Algebra, every polynomial of degree d has d complex roots. You can find these complex roots by using the POLYROOT function in SAS IML.

The programmer only wanted to output the real roots for each polynomial. No problem! You can compute the complex roots and filter out those whose imaginary parts are nonzero. This article shows how to find only the real roots of a real polynomial in SAS. The program is demonstrated by using cubic polynomials, but the logic works for higher-degree polynomials.

Find the complex roots of a polynomial

As stated earlier, the POLYROOT function in SAS IML returns the complex roots of a real polynomial of degree d. You specify the d+1 coefficients and the POLYROOT function returns the d complex roots as a d x 2 matrix where each row is a complex root. The first column is the real part of the root; the second column is the imaginary part. You must specify the coefficients in decreasing order from the degree-d mononomial down to the constant terms. You must specify a 0 for any monomial terms that do not appear in the polynomial.

For example, the following call to PROC IML computes the complex roots of the polynomial p(x) = x^3 - 3*x^2 - x + 3. Note that the vector of coefficients includes '1' for the coefficient of x3 and x. The graph of the polynomial and its roots are displayed in the graph at the top of this article.

proc iml;
/* cubic polynomial that has three real roots */
c = {1 -3 -1  3};        /* p(x) = 1*x^3 - 3*x^2 - 1*x + 3 */
roots = polyroot( c );   /* returns all complex roots */
 
Label = 'root1':'root3';
print roots[rowname=Label colname={'Re' 'Im'}];

This polynomial has three real roots, so every element in the second column is 0. In general, some polynomials will have complex roots. To get only the real roots, we can use the LOC function to select only the rows for which the second column is 0, as follows:

realIdx = loc( roots[,2] = 0 );           /* get only the real roots */
numRealRts = ncol(realIdx);               /* how many real roots? */
realRoots = rowvec( roots[realIdx, 1] );  /* put the real parts in a row vector */
print realRoots[colname=Label];
QUIT;

Although the POLYROOT function returns a matrix where each row is a complex root, it is more convenient to put the real roots in a row vector where each column is a real root.

Find the real roots of a polynomial

The previous section shows the main idea: Use POLYROOT to get the complex roots, then look at the imaginary terms and extract the real roots. However, there are a few weaknesses in the previous program:

  1. You should never perform floating point comparison with 0.
  2. The number of real roots might be less than the degree of the polynomial.
  3. It would be preferable to encapsulate the logic into a reusable function.

The following call to PROC IML implements a first version of a function that accepts a vector of coefficients and returns the real roots of the associated polynomial.

proc iml;
/* first attempt at returning the real roots of a polynomial */
start RealPolyRoot1( c, delta=1E-14 );
  roots = polyroot( c );     /* returns all complex roots */
  deg = ncol(c) - 1;         /* the degree of the polynomial */
  realRoots = j(1, deg, .);  /* there can be deg real roots */
 
  /* get only the real roots */
  realIdx = loc( abs(roots[,2]) < delta );    /* state that x + i*y is real when |y| < delta */
  numRealRts = ncol(realIdx);                 /* how many real roots? */
  realRoots[,1:numRealRts] = rowvec( roots[realIdx, 1] );  /* put the real parts in a row vector */
  return( realRoots );
finish;
 
/* test the function on two examples */
/* a cubic polynomial that has three real roots */
c = {1 -3 -1  3};        /* p(x) = 1*x^3 - 3*x^2 - 1*x + 3 */
realRoots = RealPolyRoot1( c );
Label = 'root1':'root3';
print realRoots[colname=Label];
 
/* a cubic polynomial that has one real root */
c = {1  -0.5  1  -0.5};        /* p(x) = 1*x^3 -0.5*x^2 + 1*x - 0.5 */
realRoots = RealPolyRoot1( c );
print realRoots[colname=Label];

For both examples, the function returns a three-element row vector that contains the real roots. For the second example, there is only one real root. Therefore, the vector of roots contains two missing values.

A SAS function to find the real roots of a polynomial

The function in the previous section can be improved in the following ways:

  • The function can accept a matrix of coefficients, where each row contains the coefficients of a different polynomial. In that case, the function should return a matrix of roots, where the kth row contains the roots of the kth row of coefficients.
  • The function can support coefficients that have leading zeros.

The second requirement might not be needed for your application, but it is nice to have. It is useful if some polynomials are quartic, some are cubic, some are quadratic, and so forth. The POLYROOT function requires that the leading coefficient be nonzero. Consequently, it is useful to create a helper function that strips off leading zero coefficients in a vector. You can then call the helper function before calling POLYROOT.

The following program uses the LOC function to implement the helper function. It then modifies the function in the previous section to allow a matrix of coefficients. The function loops over the rows of the coefficient matrix and returns a matrix of real roots.

proc iml;
/* Helper function: remove leading zeros or missing value from the front 
   of a vector. Examples;
   {1 0 2 3} -> {1 0 2 3}
   {0 0 1 2 3} -> {1 2 3}
*/
start RemoveLeadingZeros(c);
   nzIdx = loc(c^=0 & c^=.);           /* find all nonzero indices */
   if ncol(nzIdx) = 0 then return(0);  /* handle the "zero polynomial" */
   return( c[ , nzIdx[1]:ncol(c)] );   /* return only the elements after leading zeros */
finish;
 
/* RealPolyRoots: Return a row vector of length d that contains the real roots of a polynomial.
   The argument COEFF is a k x (d+1) matrix
   where k is the number of polynomials to examine, and d is the degree of the highest polynomial
   Classify a root as a 'real root' is the magnitude of the imaginary part is less than delta.
   Be sure to specify the coefficients in order from the highest degree down to the constant */
start RealPolyRoots( coeff, delta=1e-14 );
   numPolys = nrow(coeff);
   deg = ncol(coeff)-1;
   roots = j(numPolys, deg, .);  /* every polynomial of degree d has d complex roots */
   do i = 1 to numPolys;
      /* the leading coefficient must be nonzero */
      c = RemoveLeadingZeros( coeff[i,] );
      r = polyroot( c );   /* returns all complex roots */
      /* store only the real roots for which r[,2]=0 */
      realIdx = loc( abs(r[,2]) < delta );
      numRealRts = ncol(realIdx);
      if numRealRts > 0 then 
         roots[i, 1:numRealRts] = rowvec( r[realIdx, 1] );  /* get the real parts */
   end;
   return roots;
finish;
store module=(RemoveLeadingZeros RealPolyRoots);
 
/* matrix of coefficients for polynomials up to degree = 3 */
/*       c3  c2  c1   c0         */
coeff = {1  -3   -1   3  ,      /* 3 real roots */
         1  -0.5  1  -0.5,      /* 1 real root */
         0   1    0  -4  ,      /* quadratic polynomial with 2 real roots */ 
         0   1    0   4  ,      /* quadratic polynomial with 0 real roots */ 
         1  -0.5 -4   2  };     /* 3 real roots */
 
realRoots = RealPolyRoots( coeff );
Label = 'root1':compress('root'+char(ncol(coeff))); /* construct the names from size of coeff */
print realRoots[colname=Label];

The program calls the new RealPolyRoots function on a matrix that contains the coefficients for different polynomials: three cubic polynomials and two quadratics. The real roots are displayed for each polynomial. Notice that the quadratic polynomial x2 + 4 does not have any real roots.

Summary

A SAS programmer wanted to find the real roots of multiple polynomials by specifying the polynomial coefficients. The POLYROOT function in PROC IML can find the complex roots of a polynomial when you specify the coefficients in order from the highest degree down to the constant term. (The leading coefficient cannot be zero.) You can use the LOC function to extract the real roots by excluding the roots that have a nonzero imaginary part. This article encapsulates these steps into a user-defined function that returns only the real roots of a polynomial.

Share

About Author

Rick Wicklin

Distinguished Researcher in Computational Statistics

Rick Wicklin, PhD, is a distinguished researcher in computational statistics at SAS and is a principal developer of SAS/IML software. His areas of expertise include computational statistics, simulation, statistical graphics, and modern methods in statistical data analysis. Rick is author of the books Statistical Programming with SAS/IML Software and Simulating Data with SAS.

Leave A Reply

Back to Top