topic Re: Could you please suggest a program to find the zeros to these two polynomials in SAS/IML Software and Matrix Computations

Could you please suggest a program to find the zeros to these two polynomials

Thomas_mp — Thu, 11 Jun 2015 16:52:41 GMT

Could you please suggest a program to find the zeros to these two polynomials:

48 = ((0.3-x)*2)/((1+x)^2) solution (using excel) is -0.78717

25 = ((0.1-x)*4)/(1.x)^2 sol is -0.65292

I have a list of 4,000 of these equations to solve. and the actual polynomial have 12 terms. Perhaps if I know how to do these simple ones, I can do the more complex ones.

Thank you !!!

Re: Could you please suggest a program to find the zeros to these two polynomials

Rick_SAS — Thu, 11 Jun 2015 17:51:43 GMT

First subtract the constant terms so that you have a polynomial of the form P(x)=0. Then find the roots of the polynomials.

Do you want only the real roots, or do you want complex roots as well?

The POLYROOT function finds all complex roots of a polynomial, but you would have to represent the functions as a12*x^12 + a11*x^11 + ... + a0. See the doc or the first example of "Finding the roots of a univariate function".

To find only the real roots, use the FROOT function. Or find all and throw out the complex roots. An example of the FROOT function is given in "A simple way to find the root of a function of one variable." For the FROOT function, you can keep the functions written as the are.

Re: Could you please suggest a program to find the zeros to these two polynomials

Rick_SAS — Thu, 11 Jun 2015 17:54:22 GMT

Um, I just reread your question. Those are not polynomials, they are rational functions.

Re: Could you please suggest a program to find the zeros to these two polynomials

Thomas_mp — Fri, 12 Jun 2015 20:20:27 GMT

Thank very much Rick.

How could use this to solve several thousands functions with same "functional form"? Please, let me explain:

In your example in the link http://blogs.sas.com/content/iml/2014/02/05/find-the-root-of-a-function.html

return( exp(-x##2) - x##3 + 5#x +1 );

let me call -2,-3, 5,1 "parameters"

I have a file with 4,000 different rows with these parameters; thus each row defines a rational function; each has (at least) one root in the interval (0.05,0.30)

I tried to read the file with the 4,000 rows (I have it in a csv format) ,

infile 'H;file.csv' ; input a b c d;

The parameters in your example are (a,b,c,d), so your example is, in a generic form: ;

( exp(-x##a) - bx##c + c#x +d );

I run this, but I get the message

z = froot("Func", intervals);

ERROR: (execution) Matrix has not been set to a value.

Why?

As you can tell my knowledge of SAS is very limited..

Thank you for your help

Re: Could you please suggest a program to find the zeros to these two polynomials

PGStats — Sat, 13 Jun 2015 02:52:29 GMT

Here is the method to find the function roots with FCMP. Note however that the SOLVE function in FCMP requires a fairly good guess as initial value to find the root.

/* Equation has the form p1 = (p2-x)*p3/(p4+x)**p5 */

data parms;

eqNo + 1;

input p1-p5;

datalines;

48 0.3 2 1 2

25 0.1 4 1 2

;

/* Define FCMP routine to find function roots */

proc fcmp outlib=sasuser.fcmp.test;

function fn(p2, p3, p4, p5, x);

if fuzz(p4+x) = 0 then return (.);

else return ((p2-x)*p3/(p4+x)**p5);

endsub;

function findZero(p1, p2, p3, p4, p5, init);

array solvopts[1] initial (0);

initial = init;

return (solve("fn", solvopts, p1, p2, p3, p4, p5, .));

endsub;

run;

options cmplib=sasuser.fcmp;

data roots;

set parms;

root = findZero(of p1-p5, -0.6); /* -0.6 is the initial guess */

run;

proc print data=roots noobs; run;

Re: Could you please suggest a program to find the zeros to these two polynomials

Rick_SAS — Sat, 13 Jun 2015 09:59:26 GMT

I think this is a difficult problem in general because root-finding works best when you have knowledge of the function and the interval on which you are searching. As PG says, a good estimate for the root is important. This problem is compounded because these rational functions are not continuous and have singularities. However, in your recent response you claim that function "has (at least) one root in the interval (0.05,0.30)."

Is it always true that the function is continuous on that interval and f(0.05) and f(0.3) have different signs? If so, then you can use the FROOT function in SAS/IML and pass in that interval. The function will return the first root that it finds in that interval.

The following example reads the parameters from a SAS data set and calls the FROOT function for each row.

data parms;
eqNo + 1;
input p1-p5;
datalines;
48 0.3 2 1 2
25 0.1 4 1 2
;

proc iml;
/* Define objective function */
start fn(x) global( p );
return( (p[2]-x)*p[3]/(p[4]+x)##p[5] - p[1] );
finish;

varNames="p1":"p5";
use parms;
read all var varNames into m;
close parms;

root = j(nrow(m), 1, .);
do i = 1 to nrow(m);
   p = m[i,];           /* store i_th row into global variable */
   /* Optional: Call function that computes interval that has root */
   interval = {-0.9 2};
   root = froot("fn", interval); /* find a root in the interval [0,1] */
end;

print root;

Message was edited by: Rick Wicklin

Re: Could you please suggest a program to find the zeros to these two polynomials

Ksharp — Sat, 13 Jun 2015 12:54:02 GMT

How about using Genetic Algorithms ? which can be used under the situation that function is not differentiable or continuous .

And you define a width range of interval to search root . Unlike Rick and PG need to define a very good start point .

I love GA a hell of a lot .

Code: Program


proc iml;

start func(x);
 v=abs((0.3-x)*2/(1+x)**2 - 48);
 return (v);
finish;

id=gasetup(1,1,1234);
call gasetobj(id,0,"func");
call gasetsel(id,100,1,0.95);
call gasetcro(id,0.4,4);
call gasetmut(id,0.4);

call gainit(id,10000,{-10,10});
niter = 100;
summary = j(niter,2);
mattrib summary [c = {"Min Value", "Avg Value"} l=""];
do i = 1 to niter;
 call garegen(id);
 call gagetval(value, id);
 summary[i,1] = value[1];
 summary[i,2] = value[:];
end;
call gagetmem(mem, value, id, 1);
print mem[ l="best member:"],
   "Min Value: " value[l = ""] ;
iteration = t(1:niter);
print iteration summary;
call gaend(id);
quit;

Xia Keshan

Re: Could you please suggest a program to find the zeros to these two polynomials

Ksharp — Sat, 13 Jun 2015 13:25:42 GMT

Rick.

Not sure . If I am wrong. Your root is not right for the first function .

I get -48.16637008 not approximate 0 .

Re: Could you please suggest a program to find the zeros to these two polynomials

Rick_SAS — Sat, 13 Jun 2015 17:26:54 GMT

Thanks! I was reading all numeric variables, so I was getting the EqNo variable and using that for the first parameter. I have corrected the program to read only the parameters.

Re: Could you please suggest a program to find the zeros to these two polynomials

Thomas_mp — Sun, 14 Jun 2015 02:37:23 GMT

Thank you again !!

This was my first experience in this forum and I am very grateful.

The 3 methods work for me. I will be using Rick's because I understand a little bit better how it works. For your information, I will be using this to solve for the implied cost of capital as explained in the paper :

Toward an Implied Cost of Capital

Author(s): William R. Gebhardt, Charles M. C. Lee and Bhaskaran Swaminathan

Source: Journal of Accounting Research, Vol. 39, No. 1 (Jun., 2001), pp. 135-176

Take care

Re: Could you please suggest a program to find the zeros to these two polynomials

Ksharp — Sun, 14 Jun 2015 04:35:20 GMT

By the way, I found a better ROOT than you for the first funtion via GA. Of course for the second function .

1 OPTIONS NONOTES NOSTIMER NOSOURCE NOSYNTAXCHECK;

59 data _null_;

60 my_root=-1.254501;

61 your_root=-0.787165;

62 my_value=(0.3-my_root)*2/(1+my_root)**2 - 48;

63 your_value=(0.3-your_root)*2/(1+your_root)**2 - 48;

64 put my_root= my_value= / your_root= your_value=;

65 run;

my_root=-1.254501 my_value=0.000085964

your_root=-0.787165 your_value=-0.00020722

NOTE: DATA statement used (Total process time):

real time 0.00 seconds

cpu time 0.00 seconds

Xia Keshan

Re: Could you please suggest a program to find the zeros to these two polynomials

PGStats — Sun, 14 Jun 2015 04:41:31 GMT

it is not a better root. The function has two roots! :smileysilly: - PG

Re: Could you please suggest a program to find the zeros to these two polynomials

Ksharp — Sun, 14 Jun 2015 05:01:52 GMT

Really ? Why ? But -1.254501 is more approximate 0 . See my LOG .

Maybe it is the problem about accuracy .

Re: Could you please suggest a program to find the zeros to these two polynomials

Ksharp — Sun, 14 Jun 2015 05:16:19 GMT

Thank you . PG .

Re: Could you please suggest a program to find the zeros to these two polynomials

Rick_SAS — Sun, 14 Jun 2015 10:31:50 GMT

As I said, there is a singularity, and the function has a root to the left of the singularity and another branch to the right.

Here is a graph of the function, which has a singularity at x = -p4.:

data rational;
p1=48; p2=0.3; p3=2; p4=1; p5=2;
/* p1=25; p2=0.1; p3=4; p4=1; p5=2; */
delta = 0.05;
do x = -3 to -p4-delta by delta;
   Branch = "Neg";
   y = (p2-x)*p3/(p4+x)**p5 - p1;
   output;
end;
do x = -p4+delta to 1 by delta;
   Branch = "Pos";
   y = (p2-x)*p3/(p4+x)**p5 - p1;
   output;
end;
run;

proc sgplot data=rational;
series x=x y=y / group=Branch;
refline 0 /axis=y;
refline -1 /axis=x;
yaxis max =100;
run;

Re: Could you please suggest a program to find the zeros to these two polynomials

PGStats — Mon, 15 Jun 2015 15:33:04 GMT

An interesting alternative is to hijack the machinery of proc nlin with a single observation. The specification of starting values is very flexible in proc nlin.

/* Equation has the form p1 = (p2-x)*p3/(p4+x)**p5 */

data parms;

eqNo + 1;

input p1-p5;

datalines;

48 0.3 2 1 2

25 0.1 4 1 2

;

data startingValues;

set parms;

parameter = "x";

do Estimate = -0.95*p4 to -0.05*sign(p4) by 0.05*sign(p4);

output;

end;

keep eqNo parameter Estimate;

run;

proc nlin data=parms noprint outest=roots(where=(_TYPE_="FINAL"));

by eqNo;

parameters / pdata=startingValues;

model p1 = (p2-x)*p3/(p4+x)**p5;

run;

proc print data=roots noobs;

var EqNo _STATUS_ _SSE_ x;

format x best12.;

run;

Re: Could you please suggest a program to find the zeros to these two polynomials

billfish — Mon, 15 Jun 2015 16:09:29 GMT

The right-hand-side of the 2 examples is the ratio of 1 polynomial over another polynomial.

You write: "the actual polynomial have 12 terms". Do both polynomials (numerator and denominator) have 12 terms? Or does the equation have a grand total of 12 parameters?

Re: Could you please suggest a program to find the zeros to these two polynomials

Thomas_mp — Tue, 16 Jun 2015 03:08:00 GMT

Hello again Rick,

Thank you gain,

I thought that your program worked for my generic problem , but I am having some trouble . In my problem one of the terms is of the form:

x/(x*(1+x)##12)

For these type of problems, when x appears also multiplying in the denominator, it seems that the program cannot find the solution. For instance, in your example, below, slightly modified (I changed the first number from 48 to 3.3058, and added the x multiplying in the denominator), the solution, a zero, is x= 0.1, but the program seems to have problems finding it.

data parms;

eqNo + 1;

input p1-p5;

datalines;

3.3058 0.3 2 1 2

;

proc iml;
/* Define objective function */
start fn(x) global( p );
return( (p[2]-x)*p[3]/x*(p[4]+x)##p[5] - p[1] );
finish;

varNames="p1":"p5";
use parms;
read all var varNames into m;
close parms;

print root;

Re: Could you please suggest a program to find the zeros to these two polynomials

PGStats — Tue, 16 Jun 2015 03:43:02 GMT

Product and division have the same evaluation priority, without parentheses they are evaluated from left to right. Your equation is equivalent to

( (p[2]-x)*p[3]/x)*(p[4]+x)##p[5] - p[1]

It should read

(p[2]-x)*p[3]/(x*(p[4]+x)##p[5]) - p[1]

Re: Could you please suggest a program to find the zeros to these two polynomials

Rick_SAS — Tue, 16 Jun 2015 10:09:12 GMT

I expect that you will encounter all kinds of problems like this. For the FROOT function, the INTERVAL argument is supposed to be an interval [a,b] such that f(a) and f(b) have different signs and f is continuous on [a,b].

You see the comment that says:

/* Optional: Call function that computes interval that has root */

It means that you should use knowledge of your problem to compute an interval [a,b] for each pair of parameters. For rational functions that have k branches, you need to decide which root you want.

(If the denominator of your expression has k-1 zeros, then the rational function has k branches.)