find roots of a quadratic equation for each obs

[PamG]

I have a data set with y and x.  Is there a way to find the 2 roots of a quadratic equation for each obs? 

The equation is y - 2*x*R  + 3*x*R**2 =0.  

Sample data is

data have;
input y x;
datalines;
1 2
2 4
5 6
8 9
;
run;

[Tom]

The quadratic equation from Wikipedia:

So first convert your equation to the standard form.

Then just calculate the various parts. -b , 2a and sqrt(b**2-4ac). 

Make sure to check if 4ac > b**2 because those cases have imaginary roots.

If it has real roots then build that by combining the parts twice, once using addition and once using subtraction.

[PamG]

Thanks for your response.  I was looking for an easier way to find the roots without having to write the equations in data step.  For another project I will be having cubic polynomials.  In that case how do I save the 3 roots in 3 variables?  I saw some examples online who are using  PROC MODEL or PROC IML.  But they use some sort of starting numbers to find the roots.

[PaigeMiller]

I was looking for an easier way to find the roots without having to write the equations in data step.

Hint: there are a gazillion different articles out on the internet regarding how to perform specific mathematical and statistical calculations, many of these articles written by @Rick_SAS . Searching for these are always a good place to start.

See Rick's 2011 article: https://blogs.sas.com/content/iml/2011/08/03/finding-the-root-of-a-univariate-function.html

Hint: in your second message, you said "I was looking for an easier way to find the roots without having to write the equations in data step. For another project I will be having cubic polynomials." Please tell us all of your needs in your first message from now on.

[sbxkoenk]

• Find the root of a function by using the SAS DATA step
  By Rick Wicklin on The DO Loop August 19, 2020
  https://blogs.sas.com/content/iml/2020/08/19/find-root-function-sas-data-step.html
• A simple way to find the root of a function of one variable
  By Rick Wicklin on The DO Loop February 5, 2014
  https://blogs.sas.com/content/iml/2014/02/05/find-the-root-of-a-function.html
• Finding the root of a univariate function
  By Rick Wicklin on The DO Loop August 3, 2011
  https://blogs.sas.com/content/iml/2011/08/03/finding-the-root-of-a-univariate-function.html
• If you use SAS Viya and SAS Studio 5.2 or later AND your site licenses SAS/IML , then there's a code snippet called:
  Find Roots of Nonlinear Equation : Enables you to find the roots of a function of one variable. Finding the root (or zero) of a function enables you to solve nonlinear equations.
• Paper 3099-2019 Find a Function Root By Genetic Algorithm
  Xia Ke Shan, iFRE Inc., Beijing, China
  Kurt Bremser, Allianz Technology GmbH, Austria (Presenter)
  https://support.sas.com/resources/papers/proceedings19/3099-2019.pdf
• root finding, proc optmodel
  https://communities.sas.com/t5/SAS-Procedures/root-finding-proc-optmodel/m-p/298259
• Halley's method for finding roots 
  By Rick Wicklin on The DO Loop August 24, 2016
  https://blogs.sas.com/content/iml/2016/08/24/halleys-method-finding-roots.html

Ciao, Koen

[PamG]

Thanks for the link. I looked up and got some code.  However I am getting only one root.

**Example 3: Solve for a: P=2*e1 + a + a*a;
data test;
input p e1;
datalines;
12 1
24 2
36 3
;
run;

/* FCMP routine to find function roots */
proc fcmp outlib=work.fcmp.test;
function fn(e1, a);
    return (2*e1 + a + a*a );
endsub;
function findZero(p, e1,init);
    array solvopts[1] initial (0);
    initial = init;
    return (solve("fn", solvopts, p, e1, .)); 
endsub;
run;

options cmplib=work.fcmp.test;
         
data roots;
set test;
format root 5.2;
root = findZero(p, e1, 0.1);
run;

[sbxkoenk]

Sorry for still another link ... but I would do it this way :

• Home > Analytics > SAS/IML> find roots of function

https://communities.sas.com/t5/SAS-IML-Software-and-Matrix/find-roots-of-function/td-p/794659

A quadratic function (equation) may have one, two, or zero roots.

• If the discriminant is less than 0, then the quadratic has no real roots. (but the quadratic equation still has two complex-valued roots)
• If the discriminant is equal to zero then the quadratic has equal roots. (two roots with the same value, called a double root)
• If the discriminant is more than zero then it has 2 distinct roots.

Ciao,

Koen

[PaigeMiller]

The link I gave does not have a solution that uses PROC FCMP. So, I don't know what code you are using or what link you got it from.

[Ksharp]

Your object funtion
P=2*e1 + a + a*a;
is monotonic , so you can't expect to get two roots.

[Ksharp]

Yes. Check  @Rick_SAS  blogs. Risk has written many blog about this topic.

For your this special topic, you could use PROC NLIN ,SAS/IML .SAS/OR to get root of function.

https://blogs.sas.com/content/iml/2022/05/04/bootstrap-nonlinear-regression.html

https://blogs.sas.com/content/iml/2015/06/08/fit-circle.html

https://blogs.sas.com/content/iml/2018/08/15/optimization-nonlinear-constraints.html

Here is some example:

First is from Rick

data have;
input y x;
datalines;
12 1
24 2
36 3
;
run;

 
proc nlin data=have plots(stats=none)=(fitplot); 
   parameters R=1;                
   model y = 2*x + R + R**2;
run;

Second is from Rick

proc iml;
start ObjF(R)  global(x, y);    /* x and y are fixed vectors of data */
   F = ssq(2*x + R + R**2 - y);            /* minimize F = SSQ differences */
   return(F);
finish;
use have;  read all var {x y};  close ;  /* read observed data */
R = {10};                           /* initial guess for (x0, y0) */
opt = {0 1};                           /* options=(min, print level) */
call nlptr(rc,result,"ObjF", r, opt);  /* find optimum (x0, y0) */
print result;
quit;

Third is from @RobPratt 

proc optmodel;
   set OBS;
   num x {OBS};
   num y {OBS};
   read data have into OBS=[_N_] x y;
   var R >= 0;
   min F = sum {i in OBS} (2*x[i]+R+R^2  - y[i])^2;
   solve;
   print R;
quit;

[Ksharp]

If R was a vector:

data have;
input y x;
datalines;
12 1
24 2
36 3
;
run;



proc optmodel;
   set OBS;
   num x {OBS};
   num y {OBS};
   read data have into OBS=[_N_] x y;
   var R{OBS} >= 0;
   min F = sum {i in OBS} (2*x[i]+R[i]+R[i]^2  - y[i])^2;
   solve ;

   print F R;
quit;

[Rick_SAS · Accepted solution]

Initially you said you wanted to solve

(3*x)*R**2 - (2*x)*R + y =0, which has polynomial coefficients {3*x, -2*x, y}

But later you use 

1*R**2 + 1*R + 2*e1 = 0 

Which is it?

In any case, to get the roots of a low-degree polynomial, use the POLYROOT function in SAS IML. Create a SAS data set that has the polynomial coefficients ordered from highest degree to lowest (the constant term).  Read the coefficients into IML and loop over the coefficients.

A polynomial of degree d always has d complex roots. If you are interested only in the real roots, then find the roots for which the imaginary part is 0 (or very tiny).

The following code assumes that you want to solve the equation in your original post. It then creates a SAS data set that has ONLY the polynomial coefficients. The IML program extracts the REAL roots.

data have;
input y x;
datalines;
-1 2
-2 3
-5 6
-8 9
;

/* root of
   (3*x)*R**2 - 2*x*R  + y = 0 
   Let a = 3*x, b= -2*x, and c=y and use the quadratic equation
*/
data Coef;
keep c2 c1 c0;
set have;
c2 = 3*x;
c1 = -2*x;
c0 = y;
run;

proc iml;
delta = 1E-14;   /* small number to use as a custoff value */
use Coef;
read all var _num_ into coeff;  /* coefficients of R^2, R, constant */
close;

numPolys = nrow(coeff);
deg = ncol(coeff)-1;
roots = j(numPolys, deg, .);  /* poly of degree d has d complex roots */
do i = 1 to numPolys;
   r = polyroot( coeff[i,] );   /* returns all complex roots */
   print i r;
   /* 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;
PolyNumber = T(1:numPolys);
RootNumber = 1:deg;
print PolyNumber roots[c=RootNumber];

You can use the same IML program to find roots of higher-degree polynomials. For example, if you want to find cubic roots, use the following DATA step to create the coefficients and use the same IML code:

data Coef;
input c3 c2 c1 c0;
datalines;
1 -1 2 1
1 -1 -4 4
1 1 -3 -3
;

[PamG]

Thank you so much Rick!!!  This works perfectly.