commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Bruce A Johnson <johns...@umbc.edu>
Subject Re: [math] Numerical derivatives in Commons Math
Date Thu, 11 Aug 2011 22:47:00 GMT
I'd be quite interested in seeing Numerical Derivatives in CM.  There are some interesting
ideas about Numerical Differentiation here:

http://www.holoborodko.com/pavel/numerical-methods/

Bruce


On Aug 11, 2011, at 6:30 PM, Patrick Meyer wrote:

> I like the idea of adding this feature. What about an abstract class that implements
DifferentiableMultivariateRealFunction and provides the method for partialDerivative (). People
could then override the partialDerivative method if they have an analytic derivative.
> 
> Here's some code that I'm happy to contribute to Commons-math. It computes the derivative
by the central difference meathod and the Hessian by finite difference. I can add this to
JIRA when it's there.
> 
> /**
>     * Numerically compute gradient by the central difference method. Override this method
>     * when the analytic gradient is available.
>     *
>     *
>     * @param x
>     * @return
>     */
>    public double[] derivativeAt(double[] x){
>        int n = x.length;
>        double[] grd = new double[n];
>        double[] u = Arrays.copyOfRange(x, 0, x.length);
>        double f1 = 0.0;
>        double f2 = 0.0;
>        double stepSize = 0.0001;
> 
>        for(int i=0;i<n;i++){
>            stepSize = Math.sqrt(EPSILON)*(Math.abs(x[i])+1.0);//from SAS manual on nlp
procedure
>            u[i] = x[i] + stepSize;
>            f1 = valueAt(u);
>            u[i] = x[i] - stepSize;
>            f2 = valueAt(u);
>            grd[i] = (f1-f2)/(2.0*stepSize);
>        }
>        return grd;
>    }
> 
>    /**
>     * Numerically compute Hessian using a finite difference method. Override this
>     * method when the analytic Hessian is available.
>     *
>     * @param x
>     * @return
>     */
>    public double[][] hessianAt(double[] x){
>        int n = x.length;
>        double[][] hessian = new double[n][n];
>        double[] gradientAtXpls = null;
>        double[] gradientAtX = derivativeAt(x);
>        double xtemp = 0.0;
>        double stepSize = 0.0001;
> 
>        for(int j=0;j<n;j++){
>            stepSize = Math.sqrt(EPSILON)*(Math.abs(x[j])+1.0);//from SAS manual on nlp
procedure
>            xtemp = x[j];
>            x[j] = xtemp + stepSize;
>            double [] x_copy = Arrays.copyOfRange(x, 0, x.length);
>            gradientAtXpls = derivativeAt(x_copy);
>            x[j] = xtemp;
>            for(int i=0;i<n;i++){
>                hessian[i][j] = (gradientAtXpls[i]-gradientAtX[i])/stepSize;
>            }
>        }
>        return hessian;
>    }
> 
> 
> On 8/11/2011 5:36 PM, Luc Maisonobe wrote:
>> Le 11/08/2011 23:27, Fran Lattanzio a écrit :
>>> Hello,
>> 
>> Hi Fran,
>> 
>>> 
>>> I have a proposal for a numerical derivatives framework for Commons
>>> Math. I'd like to add the ability to take any UnivariateRealFunction
>>> and produce another function that represents it's derivative for an
>>> arbitrary order. Basically, I'm saying add a factory-like interface
>>> that looks something like the following:
>>> 
>>> public interface UniverateNumericalDeriver {
>>>  public UnivariateRealFunction derive(UnivariateRealFunction f, int derivOrder);
>>> }
>> 
>> This sound interesting. did you have a look at Commons Nabla UnivariateDifferentiator
interface and its implementations ?
>> 
>> Luc
>> 
>>> 
>>> For an initial implementation of this interface, I propose using
>>> finite differences - either central, forward, or backward. Computing
>>> the finite difference coefficients, for any derivative order and any
>>> error order, is a relatively trivial linear algebra problem. The user
>>> will simply choose an error order and difference type when setting up
>>> an FD univariate deriver - everything else will happen automagically.
>>> You can compute the FD coefficients once the user invokes the function
>>> in the interface above (might be expensive), and determine an
>>> appropriate stencil width when they call evaluate(double) on the
>>> function returned by the aformentioned method - for example, if the
>>> user has asked for the nth derivative, we simply use the nth root of
>>> the machine epsilon/double ulp for the stencil width. It would also be
>>> pretty easy to let the user control this (which might be desirable in
>>> some cases). Wikipedia has decent article on FDs of all flavors:
>>> http://en.wikipedia.org/wiki/Finite_difference
>>> 
>>> There are, of course, many other univariate numerical derivative
>>> schemes that could be added in the future - using Fourier transforms,
>>> Barak's adaptive degree polynomial method, etc. These could be added
>>> later. We could also add the ability to numerically differentiate at
>>> single point using an arbitrary or user-defined grid (rather than an
>>> automatically generated one, like above). Barak's method and Fornberg
>>> finite difference coefficients could be used in this case:
>>> http://pubs.acs.org/doi/abs/10.1021/ac00113a006
>>> http://amath.colorado.edu/faculty/fornberg/Docs/MathComp_88_FD_formulas.pdf 
>>> 
>>> It would also make sense to add vectorial and matrix-flavored versions
>>> of interface above. These interfaces would be slightly more complex,
>>> but nothing too crazy. Again, the initial implementation would be
>>> finite differences. This would also be really easy to implement, since
>>> multivariate FD coefficients are nothing more than an outer product of
>>> their univariate cousins. The Wikipedia article also has some good
>>> introductory material on multivariate FDs.
>>> 
>>> Cheers,
>>> Fran.
>>> 
>>> ---------------------------------------------------------------------
>>> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
>>> For additional commands, e-mail: dev-help@commons.apache.org
>>> 
>>> 
>> 
>> 
>> ---------------------------------------------------------------------
>> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
>> For additional commands, e-mail: dev-help@commons.apache.org
>> 
> 
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
> 
> 


---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org


Mime
View raw message