commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Patrick Meyer <meyer...@gmail.com>
Subject Re: [math] Numerical derivatives in Commons Math
Date Fri, 12 Aug 2011 11:42:15 GMT
Thanks for the information Luc. I didn't know those existed. I'm happy 
to keep the discussion at the implementation levels.

On 8/12/2011 6:23 AM, Luc Maisonobe wrote:
> Le 12/08/2011 00:30, Patrick Meyer a écrit :
>> 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.
>
> Hi Patrick,
>
> I think we need to discuss about the API we want and then about the 
> implementation. There are already other finite differences 
> implementations in some tests for both Apache Commons Math and a 
> complete package in Apache Commons Nabla. Your code adds Hessian to 
> this which is really a good thing.
>
> Thanks,
> Luc
>
>>
>> /**
>> * 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
>

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


Mime
View raw message