Hi, I just found a useful description of Jacobian construction with an example:
http://www.mathworks.de/de/help/optim/ug/writing-objective-functions.html#brkjt9_ In our case we have a vector function. Given a time, the peak intensity will be returned as I = I0.exp(-R.t). Regards, Edward On 25 August 2014 17:51, Edward d'Auvergne <[email protected]> wrote: > They should be the same ordering as the parameter vector, as each > column is a partial derivative with respect to that parameter. The > rows should be the x data or the time points. Well, I hope this is > correct, it's been many years since I coded the model-free Jacobian > matrix. Rather than using the new C code, maybe a test data set > should be calculated by hand to determine if the covariance matrix > calculation, and Jacobian calculation, is correct. > > Regards, > > Edward > > > > On 25 August 2014 17:47, Troels Emtekær Linnet <[email protected]> wrote: >> Is the order of the columns in the Jacobian matrix important? >> >> Best >> Troels >> >> 2014-08-25 17:42 GMT+02:00 Edward d'Auvergne <[email protected]>: >>> That looks correct. If you calculate: >>> >>> linalg.inv(dot(transpose(mat), mat)) >>> >>> Do you get the covariance matrix? >>> >>> Regards, >>> >>> Edward >>> >>> >>> >>> On 25 August 2014 17:35, Troels Emtekær Linnet <[email protected]> >>> wrote: >>>> Let me exemplify: >>>> >>>> There is 4 time points. >>>> >>>> df_d_i0 = - ( 2. * ( self.values - i0 / exp(r2eff * self.relax_times) >>>> ) ) / ( exp(r2eff * self.relax_times) * self.errors**2) >>>> df_d_r2eff = (2 * i0 * self.relax_times * (self.values - i0 / >>>> exp(r2eff * self.relax_times) ) ) / ( exp(r2eff * self.relax_times) * >>>> self.errors**2) >>>> >>>> Should the return then be: >>>> >>>> print df_d_i0.shape, df_d_i0 >>>> (4,) [-0.004826723918314 -0.00033019968656 0.002366308749814 >>>> -0.000232558176186] >>>> >>>> print df_d_r2eff.shape, df_d_r2eff >>>> (4,) [ 0. 2.66126225080615 -47.678483702132965 >>>> 9.371576058231405] >>>> >>>> mat = transpose(array( [df_d_i0, df_d_r2eff ]) ) >>>> print mat.shape, mat >>>> >>>> >>>> (4, 2) [[ -4.826723918313830e-03 0.000000000000000e+00] >>>> [ -3.301996865596296e-04 2.661262250806150e+00] >>>> [ 2.366308749814298e-03 -4.767848370213297e+01] >>>> [ -2.325581761857821e-04 9.371576058231405e+00]] >>>> >>>> >>>> Best >>>> Troels >>>> >>>> 2014-08-25 17:27 GMT+02:00 Edward d'Auvergne <[email protected]>: >>>>> Hi, >>>>> >>>>> I may have explained this incorrectly earlier: >>>>> >>>>> - The vector of partial derivatives with respect to the chi-squared >>>>> equation is the gradient. >>>>> - The vector of partial derivatives with respect to the exponential >>>>> function is the Jacobian. >>>>> >>>>> The equations and code for the exponential partial derivatives are the >>>>> same in both. It's just that they are used differently. Does this >>>>> help? >>>>> >>>>> Regards, >>>>> >>>>> Edward >>>>> >>>>> On 25 August 2014 17:26, Troels Emtekær Linnet <[email protected]> >>>>> wrote: >>>>>> Hi Edward. >>>>>> >>>>>> When writing the Jacobian, do you then derivative according to ( i0 * >>>>>> exp( -times * r2eff) ) >>>>>> or do you do the derivative according to chi2 function? >>>>>> >>>>>> I have a little hard time to figure out the code text. >>>>>> >>>>>> From minfx: >>>>>> @keyword func: The function which returns the value. >>>>>> @type func: func >>>>>> >>>>>> So, this is the chi2 function. >>>>>> >>>>>> @keyword dfunc: The function which returns the gradient. >>>>>> @type dfunc: func >>>>>> >>>>>> So, this must be the derivative of the chi2 function? >>>>>> >>>>>> So in essence. >>>>>> >>>>>> Does minfx expect a "dfunc" function which calculate the: >>>>>> >>>>>> one gradient chi2 value, subject to the input parameters? >>>>>> >>>>>> Or >>>>>> A jacobian matrix of the form: >>>>>> m X n matrix, where m is the number of time elements and n is number >>>>>> of parameters = 2. >>>>>> >>>>>> Best >>>>>> Troels >>>>>> >>>>>> 2014-08-25 15:52 GMT+02:00 Edward d'Auvergne <[email protected]>: >>>>>>> Hi Troels, >>>>>>> >>>>>>> Please see below: >>>>>>> >>>>>>> On 25 August 2014 13:01, Troels Emtekær Linnet <[email protected]> >>>>>>> wrote: >>>>>>>> 2014-08-25 11:19 GMT+02:00 Edward d'Auvergne <[email protected]>: >>>>>>>>> Hi Troels, >>>>>>>>> >>>>>>>>> Unfortunately you have gone ahead an implemented a solution without >>>>>>>>> first discussing or planning it. Hence the current solution has a >>>>>>>>> number of issues: >>>>>>>>> >>>>>>>>> 1) Target function replication. The solution should have reused the >>>>>>>>> C modules. The original Python code for fitting exponential curves >>>>>>>>> was converted to C code for speed >>>>>>>>> (http://gna.org/forum/forum.php?forum_id=1043). Note that two point >>>>>>>>> exponentials that decay to zero is not the only way that data can be >>>>>>>>> collected, and that is the reason for Sebastien Morin's >>>>>>>>> inversion-recovery branch (which was never completed). Anyway, the >>>>>>>>> code duplication is not acceptable. If the C module is extended with >>>>>>>>> new features, such as having the true gradient and Hessian functions, >>>>>>>>> then the Python module will then be out of sync. And vice-versa. If >>>>>>>>> a bug is found in one module and fixed, it may still be present in the >>>>>>>>> second. This is a very non-ideal situation for relax to be in, and is >>>>>>>>> the exact reason why I did not allow the cst branch to be merged back >>>>>>>>> to trunk. >>>>>>>> >>>>>>>> Hi Edward. >>>>>>>> >>>>>>>> I prefer not to make this target function dependent on C-code >>>>>>>> compilation. >>>>>>>> >>>>>>>> Compilation of code on windows can be quite a hairy thing. >>>>>>>> >>>>>>>> For example see: >>>>>>>> http://wiki.nmr-relax.com/Installation_windows_Python_x86-32_Visual_Studio_Express_for_Windows_Desktop#Install_Visual_Studio_Express_2012_for_Windows_Desktop >>>>>>>> >>>>>>>> Visua Studio Express is several hundreds of megabyte installation, for >>>>>>>> just compiling an exponential curve. ? >>>>>>>> This is way, way overkill for this situation. >>>>>>> >>>>>>> The C code compilation has been a requirement in relax since 2006. >>>>>>> This was added not only for speed, but as a framework to copy for >>>>>>> other analysis types in the future. Once a Python target function has >>>>>>> been fully optimised, for the last speed up the code can be converted >>>>>>> to C. This is the future plan for a number of the relax analyses. >>>>>>> But first the Python code is used for prototyping and for finding the >>>>>>> fastest implementation/algorithm. >>>>>>> >>>>>>> The C compilation will become an even greater requirement once I write >>>>>>> C wrapper code for QUADPACK to eliminate the last dependencies on >>>>>>> Scipy. And the C compilation framework allows for external C and >>>>>>> FORTRAN libraries to be added to the 'extern' package in the future, >>>>>>> as there are plenty of open source libraries out there with compatible >>>>>>> licences which could be very useful to use within relax. >>>>>>> >>>>>>> >>>>>>>>> 2) Scipy is now a dependency for the dispersion analysis! Why was >>>>>>>>> this not discussed? Coding a function for calculating the covariance >>>>>>>>> matrix is basic. Deriving and coding the real gradient function is >>>>>>>>> also basic. I do not understand why Scipy is now a dependency. I >>>>>>>>> have been actively trying to remove Scipy as a relax dependency and >>>>>>>>> only had a single call for numeric quadratic intergration via QUADPACK >>>>>>>>> wrappers left to remove for the frame order analysis. Now Scipy is >>>>>>>>> back :( >>>>>>>> >>>>>>>> Hi Edward. >>>>>>>> >>>>>>>> Scipy is a dependency for trying calculation with >>>>>>>> scipy.optimize.leastsq. >>>>>>>> >>>>>>>> How could it be anymore different? >>>>>>>> >>>>>>>> What you are aiming at, is to add yet another feature for estimating >>>>>>>> the errors. >>>>>>>> A third solution. >>>>>>>> >>>>>>>> What ever the third solution would come up with of dependency, would >>>>>>>> depend on the method implemented. >>>>>>>> One could also possible imagine to extend this procedures in R, Matlab >>>>>>>> or whatever. >>>>>>>> >>>>>>>> Byt they would also need to meet some dependencies. >>>>>>>> >>>>>>>> Of course the best solution would always try to make relax most >>>>>>>> independent. >>>>>>>> >>>>>>>> But if the desire is to try with scipy.optimize.leastsq, then you are >>>>>>>> bound with this dependency. >>>>>>> >>>>>>> That's why I asked if only the covariance matrix is required. Then we >>>>>>> can replace the use of scipy.optimize.leastsq() with a single function >>>>>>> for calculating the covariance matrix. >>>>>>> >>>>>>> >>>>>>>>> 3) If the covariance function was coded, then the specific analysis >>>>>>>>> API could be extended with a new covariance method and the >>>>>>>>> relax_disp.r2eff_estimate user function could have simply been called >>>>>>>>> error_estimate.covariance_matrix, or something like that. Then this >>>>>>>>> new error_estimate.covariance_matrix user function could replace the >>>>>>>>> monte_carlo user functions for all analyses, as a rough error >>>>>>>>> estimator. >>>>>>>> >>>>>>>> That would be the third possibility. >>>>>>> >>>>>>> ..., that would give the same result, save the same amount of time, >>>>>>> but would avoid the new Scipy dependency and be compatible with all >>>>>>> analysis types ;) >>>>>>> >>>>>>> >>>>>>>>> 4) For the speed of optimisation part of the new >>>>>>>>> relax_disp.r2eff_estimate user function, this is not because scipy is >>>>>>>>> faster than minfx!!! It is the choice of algorithms, the numerical >>>>>>>>> gradient estimate, etc. >>>>>>>>> (http://thread.gmane.org/gmane.science.nmr.relax.scm/22979/focus=6812). >>>>>>>> >>>>>>>> This sound good. >>>>>>>> >>>>>>>> But I can only say, that as I user I meet a "big wall of time >>>>>>>> consumption", for the error >>>>>>>> estimation of R2eff via Monte-Carlo. >>>>>>>> >>>>>>>> As a user, I needed more options to try out. >>>>>>> >>>>>>> The idea of adding the covariance matrix error estimate to relax is a >>>>>>> great idea. Despite its lower quality, it is hugely faster than Monte >>>>>>> Carlo simulations. It has been considered it before, see >>>>>>> http://thread.gmane.org/gmane.science.nmr.relax.user/602/focus=629 and >>>>>>> the discussions in that thread. But the time required for Monte Carlo >>>>>>> simulations was never an issue so the higher quality estimate remained >>>>>>> the only implementation. >>>>>>> >>>>>>> What I'm trying to do, is to direct your solution to be general and >>>>>>> reusable. I'm also thinking of other techniques at the same time, >>>>>>> Jackknife simulations for example, which could be added in the future >>>>>>> by developers with completely different interests. >>>>>>> >>>>>>> >>>>>>>>> 5) Back to Scipy. Scipy optimisation is buggy full stop. The >>>>>>>>> developers ignored my feedback back in 2003. I assumed that the >>>>>>>>> original developers had just permanently disappeared, and they really >>>>>>>>> never came back. The Scipy optimisation code did not change for many, >>>>>>>>> many years. While it looks like optimisation works, in some cases it >>>>>>>>> does fails hard, stopping in a position in the space where there is no >>>>>>>>> minimum! I added the dx.map user function to relax to understand >>>>>>>>> these Scipy rubbish results. And I created minfx to work around these >>>>>>>>> nasty hidden failures. I guess such failures are due to them not >>>>>>>>> testing the functions as part of a test suite. Maybe they have fixed >>>>>>>>> the bugs now, but I really can no longer trust Scipy optimisation. >>>>>>>>> >>>>>>>> >>>>>>>> I am sorry to hear about this. >>>>>>>> >>>>>>>> And I am totally convinced that minfx is better for minimising the >>>>>>>> dispersion models. >>>>>>>> You have proven that quite well in your papers. >>>>>>>> >>>>>>>> I do though have a hard time believing that minimisation of an >>>>>>>> exponential function should be >>>>>>>> subject to erroneous results. >>>>>>>> >>>>>>>> Anyway, this is still left to "freedom of choice" for the user. >>>>>>> >>>>>>> The error in the original Scipy optimisation code was causing quite >>>>>>> different results. The 3 algorithms, now that I look back at my >>>>>>> emails from 2003, are: >>>>>>> >>>>>>> - Nelder-Mead simplex, >>>>>>> - Levenberg-Marquardt, >>>>>>> - NCG. >>>>>>> >>>>>>> These are still all present in Scipy, though I don't know if the code >>>>>>> is different from back in 2003. The error in the Levenberg-Marquardt >>>>>>> algorithm was similar to the Modelfree4 problem, in that a lamba >>>>>>> matrix updating condition was incorrectly checked for. When the >>>>>>> gradient was positive, i.e. up hill, the matrix should update and the >>>>>>> algorithm continue to try to find a downhill step. If the conditions >>>>>>> are not correctly checked for, the algorithm thinks that the up hill >>>>>>> step means that it is at the minimum. But this is not the case, it is >>>>>>> just pointing in the wrong direction. I don't remember what the NCG >>>>>>> bug was, but that one was much more severe and the results were >>>>>>> strange. >>>>>>> >>>>>>> Failures of optimisation algorithms due to bugs can be quite random. >>>>>>> And you often don't see them, as you don't know what the true result >>>>>>> really is. But such bugs will affect exponential functions, despite >>>>>>> their simplicity. >>>>>>> >>>>>>> Regards, >>>>>>> >>>>>>> Edward _______________________________________________ relax (http://www.nmr-relax.com) This is the relax-devel mailing list [email protected] To unsubscribe from this list, get a password reminder, or change your subscription options, visit the list information page at https://mail.gna.org/listinfo/relax-devel
