diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..9993c21 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +* alias +.DS_Store diff --git a/readme.md b/readme.md index 6fceec6..05d8dd2 100644 --- a/readme.md +++ b/readme.md @@ -1,17 +1,52 @@ # Goal -This project implements a Matlab/Octave non-intrusive forward automatic differentiation method, ([wikipedia definition here](https://en.wikipedia.org/wiki/Automatic_differentiation#Forward_accumulation)) based on operator overloading. This does not provide backward mode or higher order derivatives. More specifically it implements the method refered as *sparse batched AD* in Griewank and Walther [2] Chapter 7, which allows to exploit sparsity in the intermediate Jacobians to make the computation more efficient. It enables precise and efficient computation of the Jacobian of a function. This contrasts with numerical differentiation (a.k.a finite differences) that is unprecise due to roundoff errors and that cannot exploit the sparsity of the derivatives. - -In contrast with most existing automatic differentiation Matlab toolboxes: - - * Derivatives are represented as sparse matrices, which yield to large speedups with respect to other forward mode methods when the Jacobian of the function we aim to differentiate is sparse or when intermediate accumulated Jacobian matrices are sparse (see the image denoising example). - * N dimensional arrays are supported while many Matlab automatic differentiation toolboxes only support scalars, vectors and 2D matrices - -It is likely that the speed could be improved by representing Jacobian matrices by their transpose, due to the way Matlab represents internally sparse matrices. The document [1] describes a method similar to the one implemented here and could be a very valuable source to improve the code. - -It has been tested on Matlab 2014a and Octave 4.0.0, but the example using the anonymous function @(x)eig(x) does not work on octave as octave does not call the overloaded eig function once anonymized. - -Note that backward differentation (a.k.a. gradients back-propagation in deep learning) is going to me much faster than forward differentiation when the dimension of the output is small in comparison to the dimension of the input. Forward differentiation is of interest when solving a non linear least squares for examples through Levenberg-Marquardt minimization where we want to compute the full jacobian matrix of the residuals. +This project implements a Matlab/Octave non-intrusive +forward automatic differentiation method, ([wikipedia +definition here] (https://en.wikipedia.org/wiki/ +Automatic_differentiation#Forward_accumulation)) based on +operator overloading. This does not provide backward mode or +higher order derivatives. More specifically it implements +the method refered as *sparse batched AD* in Griewank and +Walther [2] Chapter 7, which allows to exploit sparsity in +the intermediate Jacobians to make the computation more +efficient. It enables precise and efficient computation of +the Jacobian of a function. This contrasts with numerical +differentiation (a.k.a finite differences) that is unprecise +due to roundoff errors and that cannot exploit the sparsity +of the derivatives. + +In contrast with most existing automatic differentiation +Matlab toolboxes: + + * Derivatives are represented as sparse matrices, which + yield to large speedups with respect to other forward mode + methods when the Jacobian of the function we aim to + differentiate is sparse or when intermediate accumulated + Jacobian matrices are sparse (see the image denoising + example). * N dimensional arrays are supported while many + Matlab automatic differentiation toolboxes only support + scalars, vectors and 2D matrices + +It is likely that the speed could be improved by +representing Jacobian matrices by their transpose, due to +the way Matlab represents internally sparse matrices. The +document [1] describes a method similar to the one +implemented here and could be a very valuable source to +improve the code. + +It has been tested on Matlab 2014a and Octave 4.0.0, but the +example using the anonymous function @(x)eig(x) does not +work on octave as octave does not call the overloaded eig +function once anonymized. + +Note that backward differentation (a.k.a. gradients +back-propagation in deep learning) is going to me much +faster than forward differentiation when the dimension of +the output is small in comparison to the dimension of the +input. Forward differentiation is of interest when solving a +non linear least squares for examples through +Levenberg-Marquardt minimization where we want to compute +the full jacobian matrix of the residuals. # Licence @@ -65,7 +100,9 @@ more examples can be found in [./src/AutoDiffExamples.m](./src/examplesSmall.m) 0 0 0 2.0000 0 0 0 2.0000 ``` -* a simple images denoising example using a total variation (TV) regularization can be found in [./src/AutoDiffExamples.m](./src/exampleDenoise.m) +* a simple images denoising example using a total variation +(TV) regularization can be found in +[./src/AutoDiffExamples.m] (./src/exampleDenoise.m) ```c f=@(x) sum(x.^2,3); AutoDiffJacobianFiniteDiff(f,ones(2,2,2)) @@ -91,7 +128,8 @@ more examples can be found in [./src/AutoDiffExamples.m](./src/examplesSmall.m) ``` -* a simple SVM classifier training example can be found in [./src/AutoDiffExamples.m](./src/exampleSVM.m) +* a simple SVM classifier training example can be found in +[./src/AutoDiffExamples.m] (./src/exampleSVM.m) ```c function exampleSVM() % create some fake data @@ -135,21 +173,41 @@ more examples can be found in [./src/AutoDiffExamples.m](./src/examplesSmall.m) end ``` # Related projects -* [Autodiff_R2016b](https://uk.mathworks.com/matlabcentral/fileexchange/61849-autodiff_r2016b) and [Autodiff_R2015b](http://mathworks.com/matlabcentral/fileexchange/56856-autodiff) by Ultrich Reif. It uses cells to represent derivatives and uses loops instead of vectorized operations in some of the functions, which may make it too slow when using large matrices. - -* [TOMLAB/MAD](http://tomopt.com/tomlab/products/mad/). Not free. Method described in [1]. Like our code it uses operator overloading and can use sparse matrices to store directional derivatives. - -* [Automatic Differentiation for Matlab](http://www.mathworks.com/matlabcentral/fileexchange/15235-automatic-differentiation-for-matlab/) by Martin Fink. - Forward mode AD using operator overloading. Does not work with ND arrays. Not efficient for functions with sparse jacobians as it uses dense 3D arrays to store the derivatives. - -* [Automatic Differentiation with Matlab Objects](http://mathworks.com/matlabcentral/fileexchange/26807-automatic-differentiation-with-matlab-objects) by William Mcllhagga. Supports sparse jacobians but does not support ND arrays or even some matrix operations. This will fail. +* [Autodiff_R2016b] +(https://uk.mathworks.com/matlabcentral/fileexchange/61849- +autodiff_r2016b) and [Autodiff_R2015b] +(http://mathworks.com/matlabcentral/fileexchange/56856- +autodiff) by Ultrich Reif. It uses cells to represent +derivatives and uses loops instead of vectorized operations +in some of the functions, which may make it too slow when +using large matrices. + +* [TOMLAB/MAD] (http://tomopt.com/tomlab/products/mad/). Not +free. Method described in [1]. Like our code it uses +operator overloading and can use sparse matrices to store +directional derivatives. + +* [Automatic Differentiation for Matlab] +(http://www.mathworks.com/matlabcentral/fileexchange/15235- +automatic-differentiation-for-matlab/) by Martin Fink. +Forward mode AD using operator overloading. Does not work +with ND arrays. Not efficient for functions with sparse +jacobians as it uses dense 3D arrays to store the +derivatives. + +* [Automatic Differentiation with Matlab Objects] +(http://mathworks.com/matlabcentral/fileexchange/26807- +automatic-differentiation-with-matlab-objects) by William +Mcllhagga. Supports sparse jacobians but does not support ND +arrays or even some matrix operations. This will fail. ```c f=@(x) sum(x'*x) [x,dx] = autodiff(rand(5,1),f) ``` -* [madiff](https://github.com/gaika/madiff) - Reverse mode AD using operator overloading. Operators like transpose are not coded yet at the date of july 2016 . This will fail +* [madiff] (https://github.com/gaika/madiff) Reverse mode AD +using operator overloading. Operators like transpose are not +coded yet at the date of july 2016 . This will fail ```c f=@(x)(sum(x'*x)) f(rand(20,1)) @@ -158,25 +216,46 @@ more examples can be found in [./src/AutoDiffExamples.m](./src/examplesSmall.m) ``` -* [AD_deriv](https://github.com/jborggaard/AD_Deriv) by Jeff Borggaard - works only with scalars at the date of july 2016 (no vector , matrices and NDarrays) +* [AD_deriv] (https://github.com/jborggaard/AD_Deriv) by +Jeff Borggaard works only with scalars at the date of july +2016 (no vector , matrices and NDarrays) -* [Sparsegrad](https://pypi.org/project/sparsegrad/) by Marek Szymanski. Python. Automatically and efficiently calculates analytical sparse Jacobian of arbitrary numpy vector valued functions. Does not support ND arrays yet in August 2019. +* [Sparsegrad] (https://pypi.org/project/sparsegrad/) by +Marek Szymanski. Python. Automatically and efficiently +calculates analytical sparse Jacobian of arbitrary numpy +vector valued functions. Does not support ND arrays yet in +August 2019. -* [PTNobel/AutoDiff](https://github.com/PTNobel/AutoDiff) By Part Nobel. Python. Non-intrusive Forward differentiation with sparse Jacobians support. +* [PTNobel/AutoDiff] (https://github.com/PTNobel/AutoDiff) +By Part Nobel. Python. Non-intrusive Forward differentiation +with sparse Jacobians support. ## Projects that use this library -* [pde1dm](https://github.com/wgreene310/pde1dm). 1D Partial Differential Equation Solver for MATLAB and Octave. -* [NSCool_Old](https://github.com/Axect/NSCool_Old). Neutron Star Cooling simulation. -* [MatlabGP](https://github.com/noblec04/MatlabGP).A Gaussian Process and Neural Network library that is using AutoDiff for hyperparameters finding and for providing gradients for Bayesian Optimization (BO) adaptive sampling. -* [EquationFreeGit](https://github.com/uoa1184615/EquationFreeGit). Equation-Free function toolbox for Matlab/Octave. +* [pde1dm] (https://github.com/wgreene310/pde1dm). 1D +Partial Differential Equation Solver for MATLAB and Octave. + +* [NSCool_Old] (https://github.com/Axect/NSCool_Old). +Neutron Star Cooling simulation. -Please add a comment in [this issue](https://github.com/martinResearch/MatlabAutoDiff/issues/16) if you which to add you project in this listing. I am very interested in knowing what it has been used for. +* [MatlabGP] (https://github.com/noblec04/MatlabGP).A +Gaussian Process and Neural Network library that is using +AutoDiff for hyperparameters finding and for providing +gradients for Bayesian Optimization (BO) adaptive sampling. + +Please add a comment in [this issue] +(https://github.com/martinResearch/MatlabAutoDiff/issues/16) +if you which to add you project in this listing. I am very +interested in knowing what it has been used for. ## References -* [1] Forth, Shaun A. *An Efficient Overloaded Implementation of Forward Mode Automatic Differentiation in MATLAB* -ACM Trans. Math. Softw. 2006 [pdf](https://core.ac.uk/download/files/23/139791.pdf) -* [2] Andreas Griewank and Andrea Walther. *Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation*, Second Edition. 2008 +* [1] Forth, Shaun A. *An Efficient Overloaded +Implementation of Forward Mode Automatic Differentiation in +MATLAB* ACM Trans. Math. Softw. 2006 +[pdf](https://core.ac.uk/download/files/23/139791.pdf) + +* [2] Andreas Griewank and Andrea Walther. *Evaluating +Derivatives: Principles and Techniques of Algorithmic +Differentiation*, Second Edition. 2008 diff --git a/src/AutoDiff.m b/src/AutoDiff.m index 43db120..3260fda 100644 --- a/src/AutoDiff.m +++ b/src/AutoDiff.m @@ -1,45 +1,59 @@ % License FreeBSD: % -% Copyright (c) 2016 Martin de La Gorce +% Copyright (c) 2016 Martin de La Gorce, with extensions +% by Chris Noble C.2025, and Tony Roberts (AJR) June 2026. % All rights reserved. % -% Redistribution and use in source and binary forms, with or without -% modification, are permitted provided that the following conditions are met: +% Redistribution and use in source and binary forms, with or +% without modification, are permitted provided that the +% following conditions are met: % -% 1. Redistributions of source code must retain the above copyright notice, this -% list of conditions and the following disclaimer. -% 2. Redistributions in binary form must reproduce the above copyright notice, -% this list of conditions and the following disclaimer in the documentation -% and/or other materials provided with the distribution. +% 1. Redistributions of source code must retain the above +% copyright notice, this list of conditions and the +% following disclaimer. +% 2. Redistributions in binary form must reproduce the above +% copyright notice, this list of conditions and the +% following disclaimer in the documentation and/or other +% materials provided with the distribution. % -% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND -% ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED -% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE -% DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR -% ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -% (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; -% LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND -% ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -% (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND +% CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED +% WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +% OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +% INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +% (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE +% GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR +% BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +% (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +% OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +% POSSIBILITY OF SUCH DAMAGE. % -% The views and conclusions contained in the software and documentation are those -% of the authors and should not be interpreted as representing official policies, -% either expressed or implied, of the FreeBSD Project. +% The views and conclusions contained in the software and +% documentation are those of the authors and should not be +% interpreted as representing official policies, either +% expressed or implied, of the FreeBSD Project. +% +% Sometimes info obtained by command such as "help AutoDiff.svd" + classdef AutoDiff % - % This class implement a forward automatic differentation method based - % on operator overloading. This class allows precise and efficient - % computation of function Jacobians by calling AutoDiffJacobianAutoDiff + % This class implement a forward automatic differentation + % method based on operator overloading. This class allows + % precise and efficient computation of function Jacobians + % by calling AutoDiffJacobianAutoDiff % % In contrast with most AD matlab tools % - Derivatives are represented as sparse matrices % - N dimensional array are supported % - % The speed could be improved by representing jacobian matrices by - % their transposed matrix , due to the way matlab store sparse matrices + % The speed could be improved by representing jacobian + % matrices by their transposed matrix, due to the way + % matlab store sparse matrices % @@ -92,6 +106,11 @@ 'Note that vectorizing your code is likely to avoid the preallocation is likely to lead to faster execution']); end + function y = isinf(x) + % isinf() for AutoDiff x is true for infinite values + y = isinf(x.values); + end + function x = sinh(x) x.derivatives = AutoDiff.spdiag(cosh(x.values)) * x.derivatives; x.values = sinh(x.values); @@ -103,28 +122,36 @@ end function x = asinh(x) - x.derivatives = AutoDiff.spdiag(1./sqrt(1+(x.values).^2)) * x.derivatives; + x.derivatives = (1 ./sqrt(1+x.values(:).^2)).* x.derivatives; x.values = asinh(x.values); end - function x = abs(x) - x.derivatives = AutoDiff.spdiag(sign(x.values)) * x.derivatives; - x.values = abs(x.values); + function z = abs(z) + % AJR 7/2026: warning, abs() is not analytic + % AJR 7/2026: revised for both real and complex values and derivatives + sz = sign(z.values); + z.derivatives = real(sz(:)).*real(z.derivatives) ... + +imag(sz(:)).*imag(z.derivatives) ; + z.values = abs(z.values); end - function x = acosh(x) - x.derivatives = AutoDiff.spdiag(1./sqrt((x.values).^2-1)) * x.derivatives; - x.values = acosh(x.values); + function z = acosh(z) + % AJR, 7/2026: corrected sqrt for real-part negative values + s = sign(real(z.values(:))); + s( (imag(z.values)==0) & (abs(z.values+0.5)<0.5) )=+1; + z.derivatives = (s./sqrt(z.values(:).^2-1)).* z.derivatives; + z.values = acosh(z.values); end function x = atanh(x) - x.derivatives = AutoDiff.spdiag(1./(1-(x.values).^2)) * x.derivatives; + x.derivatives = AutoDiff.spdiag(1 ./(1-(x.values).^2)) * x.derivatives; x.values = atanh(x.values); end - function x = sqrt(x) - x.values = sqrt(x.values); - x.derivatives = AutoDiff.spdiag(0.5./x.values) * x.derivatives; + function z = sqrt(z) + % AJR, 7/2026: good for complex and real-negative + z.values = sqrt(z.values); + z.derivatives = ( 0.5./z.values(:) ).* z.derivatives; end function x = cos(x) @@ -144,12 +171,12 @@ end function x = acos(x) - x.derivatives = AutoDiff.spdiag(-1./sqrt(1 - x.values.^2)) * x.derivatives; + x.derivatives = AutoDiff.spdiag(-1 ./sqrt(1 - x.values.^2)) * x.derivatives; x.values = acos(x.values); end function x = asin(x) - x.derivatives = AutoDiff.spdiag(1./sqrt(1 - x.values.^2)) * x.derivatives; + x.derivatives = AutoDiff.spdiag(1 ./sqrt(1 - x.values.^2)) * x.derivatives; x.values = asin(x.values); end @@ -161,42 +188,130 @@ y = floor(x.values); end - function x = real(x) - x.derivatives = real(x.derivatives) - x.values = real(x.values); + function z = real(z) + % AJR 7/2026: warning not analytic + % and only applies when the unknown 'perturbation variable' is real + z.derivatives = real(z.derivatives);%added ; AJR, 3 Jun 2026 + z.values = real(z.values); end - function x = imag(x) - x.derivatives = imag(x.derivatives) - x.values = imag(x.values); + function z = imag(z) + % AJR 7/2026: warning not analytic + % and only applies when the unknown 'perturbation variable' is real + z.derivatives = imag(z.derivatives);%added ; AJR, 3 Jun 2026 + z.values = imag(z.values); end function x = atan(x) - x.derivatives = AutoDiff.spdiag(1./(1 + x.values.^2)) * x.derivatives; + x.derivatives = AutoDiff.spdiag(1 ./(1 + x.values.^2)) * x.derivatives; x.values = atan(x.values); end + function x = erf(x) + % erf(x) for AutoDiff x is Error function and derivatives + tmp = 2*exp(-1*x.values.^2)/sqrt(pi); + x.derivatives = AutoDiff.spdiag(tmp) * x.derivatives; + x.values = erf(x.values); + end + + function x = erfc(x) + % erfc(x) for AutoDiff x is Complementary Error function and + % derivatives. + tmp = 2*exp(-1*x.values.^2)/sqrt(pi); + x.derivatives = -1*AutoDiff.spdiag(tmp) * x.derivatives; + x.values = erfc(x.values); + end + + function x = gamma(x) + % gamma(x) for AutoDiff x is Gamma function and derivatives + tmp = gamma(x.values).*psi(x.values); + x.values = gamma(x.values); + x.derivatives = AutoDiff.spdiag(tmp) * x.derivatives; + end + function x = exp(x) x.values = exp(x.values); x.derivatives = AutoDiff.spdiag(x.values) * x.derivatives; end - function x = log(x) - tmp = 1 ./ x.values; - x.derivatives = AutoDiff.spdiag(tmp) * x.derivatives; - x.values = log(x.values); + function z = log(z) + % AJR, 7/2026: good for complex values and negative real + z.derivatives = ( 1 ./z.values(:) ).* z.derivatives; + z.values = log(z.values); end function x = tanh(x) - x.derivatives = AutoDiff.spdiag(1./(cosh(x.values).^2)) * x.derivatives; + x.derivatives = AutoDiff.spdiag(1 ./(cosh(x.values).^2)) * x.derivatives; x.values = tanh(x.values); end - function x = conj(x) - x.values = conj(x.values); - x.derivatives = conj(x.derivatives); + function z = conj(z) + % AJR 7/2026: warning not analytic + z.values = conj(z.values); + z.derivatives = conj(z.derivatives); end + function y = fft(x,varargin) + % y=fft(x,...) for AutoDiff x gives a 1D Fourier Transform + % of x and its derivatives. Optional arguments are passed + % to the usual fft. Cater for sparse derivatives as in many + % problems, although the fft-dirn may be dense, a lot of the + % variables may have no influence (e.g., Eqn-Free Patch + % scheme), that is, the 'columns' of derivatives would be + % sparse. + % + % Original C.2025 by Chris Noble, github.com/noblec04/MatlabGP + % Vectorized by AJR 7 Jun 2026. + y.values = fft(x.values,varargin{:}); + nz = size(x.derivatives,2); + nx = numel(x.values); + inz = find(~all(x.derivatives==0)); % set of non-zero columns + ninz = numel(inz); % number of non-zero columns + y.derivatives = sparse(nx,nz); + if ninz==0, y = AutoDiff(y); return, end + if 0 % original version by CN, modified + for i = inz + dxi = reshape(full(x.derivatives(:,i)),size(x.values)); + dyi = fft(dxi,varargin{:}); + y.derivatives(:,i) = dyi(:); + end%for + else % AJR, vectorized version may be roughly 10x faster + Dxi = reshape(full(x.derivatives(:,inz)),[size(x.values) ninz]); + Dxi = fft(Dxi,varargin{:}); + y.derivatives(:,inz) = reshape(Dxi,nx,ninz); + end%if 01-option + y = AutoDiff(y); + end + + function y = ifft(x,varargin) + % y=ifft(x,...) for AutoDiff x gives a 1D inverse Fourier + % Transform of x and its derivatives. Optional arguments + % are passed to the usual fft. Cater for sparse + % derivatives. + % + % Original C.2025 by Chris Noble, github.com/noblec04/MatlabGP + % Vectorized by AJR 7 Jun 2026. + y.values = ifft(x.values,varargin{:}); + nz = size(x.derivatives,2); + nx = numel(x.values); + inz = find(~all(x.derivatives==0)); % set of non-zero columns + ninz = numel(inz); % number of non-zero columns + y.derivatives = sparse(nx,nz); + if ninz==0, y = AutoDiff(y); return, end + if 0 % original version by CN, modified + for i = inz + dxi = reshape(full(x.derivatives(:,i)),size(x.values)); + dyi = ifft(dxi,varargin{:}); + y.derivatives(:,i) = dyi(:); + end%for + else % AJR, vectorized version may be roughly 10x faster + Dxi = reshape(full(x.derivatives(:,inz)),[size(x.values) ninz]); + Dxi = ifft(Dxi,varargin{:}); + y.derivatives(:,inz) = reshape(Dxi,nx,ninz); + end%if + y = AutoDiff(y); + end + function b = isreal(x) b = isreal(x.values); end @@ -262,11 +377,9 @@ x.derivatives = sparse(1:numel(r), r(:), ones(numel(r),1)) * x.derivatives; end - function x = ctranspose(x) - x = transpose(x); - if ~isreal(x) - x = conj(x); - end + function z = ctranspose(z) + z = transpose(z); + if ~isreal(z), z = conj(z); end end function D = spdiags(B, d, m, n) @@ -303,8 +416,10 @@ end function x = diff(x, n, dim) - if n ~= 1 - error('not yet coded') + if nargin<3, dim=min(find(size(x)>1)); end %AJR, 31 May 2026 + if nargin<2, n=1; end %AJR, 28 May 2026 + if n>1 %AJR, 28 May 2026, recursion probably inefficient + x = diff(x,n-1,dim); end if issparse(x.values) warning('AutoDiff:Inefficient', 'this implementation is quite inefficent') @@ -355,8 +470,20 @@ end end - function z = sign(x) - z = sign(x.values); + function z = sign(z) + % AJR 7/2026: revised, and also for complex arguments + if isreal(z) + [m,n] = size(z.derivatives); + z.derivatives = sparse(m,n); + z.values = sign(z.values); + else % sign() for complex values + r = abs(z.values); + z.derivatives = z.derivatives./r(:) ... + -z.values(:)./r(:).^3.*( ... + real(z.values(:)).*real(z.derivatives) ... + +imag(z.values(:)).*imag(z.derivatives) ); + z.values = sign(z.values); + end end function x = subsindex(x) @@ -517,26 +644,66 @@ function x = mpower(x, n) - if numel(x) == 1 - x = x.^n; - else - if n == 1 - return - elseif n > 1 - x = mtimes(x^(n - 1), x); - else - error('not coded yet') - end + if numel(x) == 1, x = x.^n; return; end + % AJR: should include negative powers and fractional + assert(round(n)==n ... + ,'AutoDiff: not yet coded fractional powers') + assert(n>=0 ... + ,'AutoDiff: not yet coded negative powers') + % AJR, 7/2026: added following case of zero power + if n == 0, x.values=eye(size(x.values)); + [k,l] = size(x.derivatives); + x.derivatives = sparse(k,l); + elseif n == 1, return + elseif n == 2, x = mtimes(x,x); + elseif n > 2 %AJR, 7/2026: log(n) multiplies version + m = floor(n/2); + if 2*m == n, y=[]; else y=x; end + x = mpower(mtimes(x,x),m); + if ~isempty(y), x = mtimes(x,y); end + else error('AutoDiff: should not occur') end end - function x = inv(x) - x.values = inv(x.values); + function z = inv(z) + z.values = inv(z.values); + M1 = kron(speye(size(z.values, 2)), z.values); + %AJR, 7/2026, needed .' instead of ' in following + M2 = kron(z.values.', speye(size(z.values, 1))); + z.derivatives = -M2 * M1 * z.derivatives; + end + + function x = pinv(x,tol) + % pinv(x) for AutoDiff x is Moore-Penrose pseudoinverse of + % x, and derivatives. By Chris Noble C.2025 + if nargin==1 + tol=0; + end + x.values = pinv(x.values,tol); M1 = kron(speye(size(x.values, 2)), x.values); M2 = kron(x.values', speye(size(x.values, 1))); x.derivatives = -M2 * M1 * x.derivatives; end + function [y,flag] = chol(x,shape) + % chol(x,shape) for AutoDiff x is Cholesky factorization of + % x, and derivatives. By Chris Noble C.2025 + if nargin<2, shape='lower'; end + [L,flag] = chol(x.values,shape); + if flag~=0, y=L; return, end + for i = 1:size(x.derivatives,2) + A = L\reshape(x.derivatives(:,i),size(x))/L'; + n = size(A,1); + A(1:(n+1):end) = 0.5*A(1:(n+1):end); + U = tril(0*A+1); + A(U~=1) = 0; + A = L*A; + y.derivatives(:,i) = A(:); + end + y.values = L; + y = AutoDiff(y); + end + function z = mldivide(x, y) if isa(y, 'AutoDiff') if isa(x, 'AutoDiff') @@ -656,7 +823,9 @@ [i,j,k,l]=ndgrid(1:size_x(1),1:size_y(1),1:size_y(2),1:s); j2 = j+(k-1)*size_y(1)+(l-1)*size_y(1)*size_y(2); i2 = i+(k-1)*size_x(1)+(l-1)*size_x(1)*size_y(2); - v=repmat(reshape(x_values,size_x(1),size_x(2),prod(size_x(3:end))),1,size_y(2),1); + v = repmat( reshape( x_values ... + ,size_x(1),size_x(2),prod(size_x(3:end)) ) ... + ,1,size_y(2),1 ); Mx=sparse(i2(:),j2(:),v(:)); end if isa(x, 'AutoDiff') @@ -679,30 +848,118 @@ if (numel(y) == 1) z = x ./ y; return; - else - error('not yet coded') + else error('not yet coded') end end - - function x = norm(x, p) - if nargin == 1 - p = 2; - end - - if isvector(x) - x = sum(abs(x.^p)).^(1 / p); - elseif ismatrix(x) - [~, d, ~] = svd(x); - x = max(d); - else - error('not sure what matlab does in this case'); + function z = norm(z, p) + % AJR 7/2026: warning, norm() is not analytic + if nargin == 1, p = 2; end + if isvector(z) + z = sum(abs(z).^p).^(1/p); + elseif ismatrix(z) + [~, d, ~] = svd(z); + z = max(d); + else error('not sure what matlab does in this case'); end end - function [U, S, V] = svd(x) - error('not coded yet, could look at the eig implementation') - end + function [Uad, Sad, Vad] = svd(A,varargin) + % Economy-size SVD A=USV' for AutoDiff matrix A, mxn real or + % complex: gives mxk U, kxk S, nxk V and its derivatives, + % for rank k=number of non-zero singular values + % [k<=min(m,n)]. Assumes there are no coincident singular + % values (as otherwise F divides by zero). Only compute the + % economy rank k decomposition because for (multiple) + % zero-rows/columns of S the corresponding columns of U&V + % are not unique, and so their derivative is meaningless. + % + % AJR 22/6/2026, much revised from Chris Noble's MatlabGP, + % and see "Differentiating the SVD", James Townsend, 2016 + % + % AutoDiffsvdRankThresh (global) sets threshold on singular + % values for determining the rank of the SVD and its + % derivatives: singular values such that s_i/s_1 < threshold + % are zeroed, and rank reduced accordingly. Reason? + % derivatives of U&V involve division by singular values so + % are very sensitive to very small ones. + if nargin>1, warning(['AutoDiff.svd: ' ... + '2nd+ argument ignored, always economy sized']), end + global AutoDiffsvdRankThresh AutoDiffsvdCmplxWarn + if ~( exist('AutoDiffsvdRankThresh') ... + && ~isempty(AutoDiffsvdRankThresh) ) + AutoDiffsvdRankThresh = 1e-8; % guess useful default + end%if exist + if ~isreal(A.values) % warn about complex derivatives + if ~( exist('AutoDiffsvdCmplxWarn') && ... + strcmp(class(AutoDiffsvdCmplxWarn),'datetime') ) + AutoDiffsvdCmplxWarn = datetime('yesterday'); + end%if ~exist + tNow = datetime('now'); + if hours(tNow-AutoDiffsvdCmplxWarn)>1 % warn every one hour + warning(['derivatives of SVD are generally NOT' ... + ' analytic functions of complex matrices']) + AutoDiffsvdCmplxWarn = tNow; + end%if hours>1 + end%if ~isreal + [m,n] = size(A.values); + [U,S,V] = svd(A.values); + s = diag(S); + k = max(find(s>AutoDiffsvdRankThresh*s(1))); + s = s(1:k); % kx1 col.vector of singular values + U = U(:,1:k); % mxk as in T2016 + S = diag(s); % kxk matrix + V = V(:,1:k); % nxk as in T2016 + % 'rotate' each column of U&V so largest-V is real-positive + % Means computed U&V is almost-always unique. + [~,iVM]=max(abs(V)); + c=nan(1,k); + for j=1:k, c(j)=V(iVM(j),j)/abs(V(iVM(j),j)); end + U=U./c; V=V./c; % AJR: modified 15 Jul 26 to avoid conj + + N = size(A.derivatives,2); + dU = sparse(m*k,N); + dS = sparse(k*k,N); + dV = sparse(n*k,N); + % cater for sparse columns in A.derivatives + [~,jN] = find(A.derivatives); % find non-zero columns + jN = unique(jN); + N = length(jN); % redefined to number of non-zero cols + % use fast pagemtimes if available + if exist('pagemtimes')==5, Mx = @pagemtimes; + else Mx = @AutoDiff.pageMmult; + end%if exist + T = 'ctranspose'; O = 'none'; + dA = reshape(full(A.derivatives(:,jN)),m,n,N); + dP = Mx(U,T,Mx(dA,V),O); % U'*dA*V k.k.N + iDiag = find(speye(k)); + dP2D = reshape(dP,k^2,N); + dS(iDiag,jN) = real(dP2D(iDiag,:)); % diag(real(diag(dP))); + dP = dP-reshape(full(dS(:,jN)),k,k,N); + Foff=(1-eye(k))./(s'.^2-s.^2+eye(k)); % k.k + symM=@(X) X+conj(permute(X,[2 1 3])); + UdUoff=Foff.*symM(dP.*s'); % Foff.*(dP*S+S*dP'), k.k.N + VdVoff=Foff.*symM( s.*dP); % Foff.*(S*dP+dP'*S), k.k.N + dVi = Mx(V, VdVoff ); % V*( VdVoff ), n.k.N + dP2D = reshape(dP,k^2,N); + UdUoff = reshape(UdUoff,k^2,N); % UdUoff+diag(diag(dP)./s) + UdUoff(iDiag,:) = UdUoff(iDiag,:)+dP2D(iDiag,:)./s; + dUi = Mx(U, reshape(UdUoff,k,k,N) ); % U*( UdUoff ), m.k.N + if m>k, dUi=dUi+Mx(Mx(eye(m)-U*U',O,dA,O),V./s'); end + if n>k, dVi=dVi+Mx(Mx(eye(n)-V*V',O,dA,T),U./s'); end + c = nan(1,k,N); + for j=1:k + c(1,j,:) = (-1i./abs(V(iVM(j),j))).*imag(dVi(iVM(j),j,:)); + end%for j + dV(:,jN) = reshape( dVi+V.*c ,n*k,N); + dU(:,jN) = reshape( dUi+U.*c ,m*k,N); + % form results as AutoDiff entities + Uad = AutoDiff(U,dU); + Sad = AutoDiff(S,dS); + Vad = AutoDiff(V,dV); + if nargout == 1, Uad = Sad; end + end% svd() + function n = numel(x) n = numel(x.values); @@ -751,7 +1008,7 @@ if isa(x, 'AutoDiff') x = repmat_as(x, y); y = repmat_as(y, x); - x.derivatives = AutoDiff.spdiag(1./y.values) * x.derivatives - AutoDiff.spdiag(x.values./y.values.^2) * y.derivatives; + x.derivatives = AutoDiff.spdiag(1 ./y.values) * x.derivatives - AutoDiff.spdiag(x.values./y.values.^2) * y.derivatives; x.values = x.values ./ y.values; else @@ -764,7 +1021,7 @@ else x = repmat_as(x, y); y = repmat_as(y, x); - x.derivatives = AutoDiff.spdiag(1./y) * x.derivatives; + x.derivatives = AutoDiff.spdiag(1 ./y) * x.derivatives; x.values = x.values ./ y; end end @@ -792,7 +1049,7 @@ end - function x = cumsum(x, varargin) + function x = cumsum(x, varargin) val = cumsum(x.values, varargin{:}); if isvector(x.values) @@ -870,7 +1127,7 @@ end %y.derivatives(tmp,:)=x.derivatives; % slow for some - %reasons for large sparse matrices + %reasons for large sparse matrices n = numel(y.values); m = numel(x); @@ -1032,78 +1289,89 @@ end - function [V, D] = eig(C) - % Compute the eigen vector eigen values and there derivative with respect - % to each element of the input matrix. The function might be undifferentiable - % if the mutiplicity of an eigen value is more than one. - % It may no work if C is not symmetric (need to check if the formulas are still valid) - if any(any(C.values' - C.values) > eps) - error('not yet verified for non symetric matrices') - end - - n = size(C, 1); - [V, D] = eig(C.values); - lambda = diag(D); - % C.values*V==V*D - % k=1 - % C.values*V(:,k)=lambda(k)*V(:,k) - % - - l = 0; - - dV_dC = zeros(n, n, n^2); - dD_dC = zeros(n, n, n^2); - - dlambda = zeros(size(C, 1), n^2); - for j = 1:n - for i = 1:n - l = l + 1; - - Ap = sparse(i, j, 1, size(C, 1), size(C, 1)); - - - for k = 1:size(C, 1) - %dlambda(k,l)=V(i,k)*V(j,k) - dlambda(k, l) = V(:, k)' * Ap * V(:, k); - - - % B=[C-lambda(k)*eye(n,n);V(:,k)']; - % dV_dC(:,k,l)=(B'*B)^-1*B'*[dlambda(k)*V(:,k)-Ap*V(:,k);0]; - dV_dC(:, k, l) = [C.values - lambda(k) * eye(n, n); V(:, k)'] \ [dlambda(k, l) * V(:, k) - Ap * V(:, k); 0]; - % [C-lambda(k)*eye(3,3)]*dV_dC(:,k,l)+-dlambda(k)*V - - - %n=size(C.values,1); - % k=1; - % - % (C.values-lambda(k)*eye(n))*V(:,k) - - % Ap=sparse(i,j,1,size(C,1),size(C,1)); - % (Ap-dlambda(k,l)*eye(n))*V(:,k)+(C.values-lambda(k)*eye(n))*dV_dC(:,k,l) - % dV_dC(:,k,l)'*V(:,k) - % V(:,k)'*(Ap-dlambda(k,l)*eye(n))*V(:,k)+V(:,k)'*(C.values-lambda(k)*eye(n))* dV_dC(:,k,l) - - % V(:,k)'*(C.values-lambda(k)*eye(n)) - end - dD_dC(:, :, l) = diag(dlambda(:, l)); - end - end - - - if nargout == 1 - V = AutoDiff(lambda, dlambda); - else - - D = AutoDiff(D, reshape(dD_dC, numel(D), [])*C.derivatives); - V = AutoDiff(V, reshape(dV_dC, numel(D), [])*C.derivatives); - end - end - - - function x = transpose(x) - M = AutoDiff.transposeDiff(size(x)); - x.derivatives = M * (x.derivatives); - x.values = x.values'; + function [V,D] = eig(Z,varargin) + % For matrix Z, find eigenvectors V, eigenvalues D, and + % their AutoDiff differentials. V and D are analytic + % functions of perturbations to Z, as seen in their + % differentials. + % Input: Z = square matrix, real or complex, symmetric or + % not, possibly sparse, class double or AutoDiff + % Outputs: V = AutoDiff nxn matrix of eigenvectors; each + % column of V scaled so the largest magnitude + % component is precisely one (a way to ensure + % analyticity). + % D = AutoDiff nxn sparse diagonal matrix of + % eigenvalues correspoonding to V, sorted in + % increasing order of real-part. + % AJR, 20 Jul 2026 + if nargin>1, warning(['AutoDiff.eig: ' ... + '2nd+ argument(s) ignored as yet']), end + [n,n2] = size(Z); + assert(n2==n,'AutoDiff.eig: matrix not square') + n2=n*n; % use for number of elements in Z + if strcmp(class(Z),'AutoDiff') + dZ = Z.derivatives; + Z = Z.values; + else% Z is double, so set dZ sparse identity + % maybe explore assuming AD w.r.t only non-zero entries + ks = 1:n2;%find(Z); + dZ = sparse(ks,ks,1,n2,n2); + end%if strcmp + % to get e-vecs, eig requires full matrix. Perhaps later + % code eigs version for finding subset of e-vecs. + [V,D] = eig(full(Z)); + % always need to sort in a consistent order, and let's sparsify D + [d,j] = sort(diag(D),'ComparisonMethod','real'); + D = spdiags(d,0,n,n); V = V(:,j); + % ensure largest element is pure real +1 + [~,iM]=max(abs(V)); + %?? for j=1:n, V(:,j) = V(:,j)/sign(V(iM(j),j)); end + for j=1:n, V(:,j) = V(:,j)/V(iM(j),j); end + % check it + for j=1:n, assert( (abs(imag(V(iM(j),j)))<1e-14)&(real(V(iM(j),j))>0) ), end + % optionally test V,D still eigensoln + %assert(norm(Z*V-V*D,'fro')<1e-12); + + % proceed to find the eigen-differentials + dV = sparse(n2,size(dZ,2)); + dD = sparse(n2,size(dZ,2)); + % loop over all e-value & e-vector pairs + for j=1:n + Vj = V(:,j); % n-D e-vector + Dj = D(j,j); % scalar e-value + % solve (Z-Dj*I)dv-Vj*dD = -dZ*Vj for all dZ_k + % adjoin that perturbations to max-V are zero as this + % leads to eigenvectors being analytic functions of the + % parameters + Zjx = [ Z-speye(n)*Dj -Vj + sparse(1,iM(j),1,1,n+1) ]; + % vectorize product dZ_k*Vj by recoding dZ(:,ks)*Vj + % Since dZ(:,k) is stored column-wise, pre-multiply dZ + % by sparse([ diag(V_1) diag(V_2) ... diag(V_n); 0s ]) + r = rem(0:n2-1,n)+1; i = ceil((1:n2)/n); + rhs = -sparse(r,1:n2,Vj(i),n+1,n2)*dZ; + dVj = Zjx\rhs; % solve simultaneously for all dZ columns + dDj = dVj(n+1,:); % e-val diffs for all k + dVj = dVj(1:n,:); % e-vecs can add arb multiple of V_j + if 0, % omit, this enforcing normality of V+dV wrecks analyticity + dVj = dVj - Vj.*real( Vj'*dVj ); + normErr = real( Vj'*dVj ); + assert(all(abs(normErr)<1e-12)) + end%if 0/1 + % store the differentials + dV((j-1)*n+(1:n),:) = dVj; + dD((j-1)*n+j,:)= dDj; + end%for j + V = AutoDiff(V,dV); + D = AutoDiff(D,dD); + if nargout==1, V=diag(D); D=[]; end + end%function eig + + + function z = transpose(z) + M = AutoDiff.transposeDiff(size(z)); + z.derivatives = M * (z.derivatives); + z.values = z.values.'; %AJR, 7/2026: needed .' instead of ' end function x = permute(x, l) @@ -1127,10 +1395,11 @@ end - function y = det(x) - y.values = det(x.values); - y.derivatives = reshape(det(x.values).*inv(x.values)',1,[]) * x.derivatives; - y = AutoDiff(y.values, y.derivatives); + function z = det(z) + % AJR, 7/2026: complex matrices x need dot-transpose + detz = det(z.values); + z.derivatives = reshape(detz.*inv(z.values).',1,[]) * z.derivatives; + z.values = detz; end function y = vertcat(varargin) @@ -1142,16 +1411,29 @@ assert(length(varargin) == k+2); assert(strcmp(varargin{k + 1}, 'like')); x.values = ones(varargin{1:k}); - x.derivatives = zeros(numel(x.values), size(varargin{k + 2}.derivatives, 2)); + x.derivatives = sparse(numel(x.values), size(varargin{k + 2}.derivatives, 2));%change zeros to sparse zero, AJR 3/6/26 x = AutoDiff(x); end function x = zeros(varargin) + % AJR: to get invoked need to use 'like'? k = find(cellfun(@isnumeric, varargin), 1, 'last'); assert(length(varargin) == k+2); assert(strcmp(varargin{k + 1}, 'like')); x.values = zeros(varargin{1:k}); - x.derivatives = zeros(numel(x.values), size(varargin{k + 2}.derivatives, 2)); + x.derivatives = sparse(numel(x.values), size(varargin{k + 2}.derivatives, 2));%changed zeros to sparse, AJR 3/6/26 + x = AutoDiff(x); + end + + function x = nan(varargin) + % nan(...) for AutoDiff-likeness gives nan-matrix according + % to varargin, and sparse-zero derivatives. AJR, 9 Jun 2026 + k = find(cellfun(@isnumeric, varargin), 1, 'last'); + assert(length(varargin) == k+2); + assert(strcmp(varargin{k + 1}, 'like')); + x.values = nan(varargin{1:k}); + x.derivatives = sparse(numel(x.values) ... + ,size(varargin{k + 2}.derivatives, 2)); x = AutoDiff(x); end @@ -1163,6 +1445,39 @@ methods (Static) + + function C = pageMmult(A,At,B,Bt) + % Computes matrix product A*B for ensemble of matrices. + % Arguments are either two (A,B), or four (A,At,B,Bt), where + % A & B are 3-D arrays. For each l computes + % A(:,:,l)*B(:,:,l) where A&B may be transposed according to + % At & Bt: 'none' (default), 'ctranspose' or 'transpose' + % AJR, 25 Jun 2026 + if nargin==2, B=At; At='none'; Bt='none'; + else assert(nargin==4 ... + ,'AutoDiff.pageMmult: must have 2 or 4 arguments') + end%if + [a1,a2,a3]=size(A); [b1,b2,b3]=size(B); + k = max(a3,b3); + if a3>1, a=@(l) l; else a=@(l) 1; end + if b3>1, b=@(l) l; else b=@(l) 1; end + switch At + case 'none', opA=@(X) X; c1=a1; + case 'ctranspose', opA=@ctranspose; c1=a2; + case 'transpose', opA=@transpose; c1=a2; + otherwise error('AutoDiff.pageMmult: unknown A-transpose') + end%switch At + switch Bt + case 'none', opB=@(X) X; c2=b2; + case 'ctranspose', opB=@ctranspose; c2=b1; + case 'transpose', opB=@transpose; c2=b1; + otherwise error('AutoDiff.pageMmult: unknown B-transpose') + end%switch Bt + C=nan(c1,c2,k); + for l=1:k, C(:,:,l)=opA(A(:,:,a(l)))*opB(B(:,:,b(l))); end + end%function pageMmult + + function M = spDiagFromVec(v) M = sparse((1:numel(v)), (1:numel(v)), v(:)); end diff --git a/src/AutoDiffJacobianAutoDiff.m b/src/AutoDiffJacobianAutoDiff.m index 8b0c4e0..6f7a120 100644 --- a/src/AutoDiffJacobianAutoDiff.m +++ b/src/AutoDiffJacobianAutoDiff.m @@ -1,7 +1,7 @@ function [J, f] = AutoDiffJacobianAutoDiff(func, x, range) -% AutoDiffJacobianAutoDiff returns the jacobian of function evaluated at x using Automatic -% Differentiation +% AutoDiffJacobianAutoDiff returns the jacobian of function +% evaluated at x using Automatic Differentiation % % % J = AutoDiffJacobianAutoDiff(func,x) @@ -15,12 +15,14 @@ % J : jacobian matrix of function f evaluated at x % % Description: -% returns the jacobian matrix of function f evaluated at x -% if a range is specified , only the columns corresponding to the range -% indices are computed. The computation is donc using Automatic -% Differentiation (Not to be confounded with Finite Differences method) -% This is numerically more precise and, if the Jacobian is large and -% sparse this might be much faster (see examples) +% returns the jacobian matrix of function f evaluated at +% x if a range is specified , only the columns +% corresponding to the range indices are computed. The +% computation is donc using Automatic Differentiation +% (Not to be confounded with Finite Differences method) +% This is numerically more precise and, if the Jacobian +% is large and sparse this might be much faster (see +% examples) % % Examples: % @@ -55,7 +57,8 @@ end sx = size(x); nr = numel(range); - xAD = AutoDiff(x, sparse(range, (1:nr), ones(1, nr), sx(1) * sx(2), nr)); + % AJR 11 Jun 2026: changed to prod() to cater for >2-D arrays x + xAD = AutoDiff(x, sparse(range, (1:nr), ones(1, nr), prod(sx), nr)); end try @@ -63,11 +66,13 @@ f = getvalue(fAD); J = getderivs(fAD); catch exception - warning('failed while calling the function with the AutoDiff instance, trying to call it with the plain data instead to chech that works') + warning('failed while calling the function with the AutoDiff instance, trying to call it with the plain data instead to check that works') func(x); warning('It seems like to original function is ok with the plain data , The class AutoDiff needs debugging') % rethrow(exception) - [~] = func(xAD); % better to call again the function instead of using exception as it makes it possible to use matlab's stop-if -error debugging functionality + [~] = func(xAD); % better to call again the function instead of using +% exception as it makes it possible to use matlab's stop-if +% -error debugging functionality end diff --git a/src/AutoDiffJacobianFiniteDiff.m b/src/AutoDiffJacobianFiniteDiff.m index 0d64665..0d958b8 100644 --- a/src/AutoDiffJacobianFiniteDiff.m +++ b/src/AutoDiffJacobianFiniteDiff.m @@ -1,7 +1,7 @@ function [J, f] = AutoDiffJacobianFiniteDiff(func, x, range, epsilons, centered) -% AutoDiffJacobianFiniteDiff returns the jacobian of function f evaluated at x using Finite -% Differences. +% AutoDiffJacobianFiniteDiff returns the jacobian of +% function func evaluated at x using Finite Differences. % % % J = AutoDiffJacobianFiniteDiff(func,x) @@ -12,16 +12,19 @@ % x : location where the jacobian is evaluated % range : (optional) vector of integer indices in [1,numel(x)] % Outputs: -% J : jacobian matrix of function f evaluated at x +% J : jacobian matrix of function func evaluated at x +% f : value of function func evaluated at x % % Description: -% returns the jacobian matrix of function f evaluated at x -% if a range is specified , only the columns corresponding to the range -% indices are computed. The computation is donc using Finite -% Differences. This might be unprecise due to rounoff errors. -% This is provided as a way to check the validity of the Automatic differenciation -% method. Considere using AutoDiffJacobianAutoDiff instead to get faster -% and more accurate derivatives. +% returns the jacobian matrix of function f evaluated at +% x if a range is specified, only the columns +% corresponding to the range indices are computed. The +% computation is done using Finite Differences. This +% might be imprecise due to rounoff errors. This is +% provided as a way to check the validity of the +% Automatic differentiation method. Consider using +% AutoDiffJacobianAutoDiff instead to get faster and +% more accurate derivatives. % % Documentation created by Martin de La Gorce diff --git a/src/AutoDiff_tests.m b/src/AutoDiff_tests.m index 1f76dff..4e4fff0 100644 --- a/src/AutoDiff_tests.m +++ b/src/AutoDiff_tests.m @@ -1,236 +1,269 @@ isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0; -if isOctave - rand ("seed", 1) -else - rng(1) -end +% AJR: specify random seed to random or toreproduce +if 1, ADseed = floor(1e5*rand(1)); +else ADseed = 1; end +if isOctave, rand ("seed", ADseed), else rng(ADseed), end if ~exist('pagemtimes') addpath("./backports") end -% to be fixed -% f = @(x) diff(x, 1, 3); -% CheckAutoDiffJacobian(f, randn(2,2), 1e-9); +% AJR, fixed June 2026 +f = @(x) diff(x, 1, 3); +CheckAutoDiffJacobian(f, randn(2,3,4), 1e-8); +f = @(x) diff(x, 1, 2); +CheckAutoDiffJacobian(f, randn(3)+1i*rand(3), 1e-8); + +% AJR, added July 2026 +f = @(x) abs(x); +CheckAutoDiffJacobian(f, randn(2)+i*randn(2), 1e-8); + +% AJR, added July 2026 +f = @(x) sign(x); +z = randn(2)+i*randn(2); z=z+0.2*sign(z); % avoid singularity at 0 +CheckAutoDiffJacobian(f, z, 1e-8); f = @(x) real(x); -CheckAutoDiffJacobian(f, randn(2,2)+i*randn(2,2), 1e-9); +CheckAutoDiffJacobian(f, randn(2)+i*randn(2), 1e-8); f = @(x) imag(x); -CheckAutoDiffJacobian(f, randn(2,2)+i*randn(2,2), 1e-9); +CheckAutoDiffJacobian(f, randn(2)+i*randn(2), 1e-8); f = @(x) i*x; -CheckAutoDiffJacobian(f, randn(2,2)+i*randn(2,2), 1e-9); +CheckAutoDiffJacobian(f, randn(2)+i*randn(2), 1e-7); f = @(x) real(i*x); -CheckAutoDiffJacobian(f, randn(2,2)+i*randn(2,2), 1e-9); +CheckAutoDiffJacobian(f, randn(2)+i*randn(2), 1e-7); f = @(x) sum(x,3); -CheckAutoDiffJacobian(f, randn(2,2), 1e-9); +CheckAutoDiffJacobian(f, randn(2,3,5), 1e-8); +CheckAutoDiffJacobian(f, randn(2,3,5)+1i*randn(2,3,5), 1e-8); -f = @(x) mean(x,3); -CheckAutoDiffJacobian(f, randn(2,2), 1e-9); +f = @(x) mean(x,2); +CheckAutoDiffJacobian(f, randn(2,3), 1e-8); +CheckAutoDiffJacobian(f, randn(2,3)+1i*randn(2,3), 1e-8); f = @(x) cat(2, [], x); -CheckAutoDiffJacobian(f, randn(3,3), 1e-9); +CheckAutoDiffJacobian(f, randn(3,3), 1e-8); f = @(x) cat(3, [], x); -CheckAutoDiffJacobian(f, randn(3,3), 1e-9); +CheckAutoDiffJacobian(f, randn(3,3), 1e-8); f = @(x) cat(3, [],[], x); -CheckAutoDiffJacobian(f, randn(3,3), 1e-9); +CheckAutoDiffJacobian(f, randn(3,3), 1e-8); f = @(x) cat(3, x, []); -CheckAutoDiffJacobian(f, randn(3,3), 1e-9); +CheckAutoDiffJacobian(f, randn(3,3), 1e-8); f = @(x) cat(4, [], x); -CheckAutoDiffJacobian(f, randn(3,3), 1e-9); +CheckAutoDiffJacobian(f, randn(3,3), 1e-8); f = @(x) cat(4, x,[]); -CheckAutoDiffJacobian(f, randn(3,3), 1e-9); +CheckAutoDiffJacobian(f, randn(3,3), 1e-8); x = randn(3, 2, 7); f = @(y) pagemtimes(x,y); -CheckAutoDiffJacobian(f, randn(2, 5, 1), 1e-9); +CheckAutoDiffJacobian(f, randn(2, 5, 1), 1e-8); -x = randn(3, 2, 1); +x = randn(3, 2, 1)+1i*randn(3, 2, 1); f = @(y) pagemtimes(x,y); -CheckAutoDiffJacobian(f, randn(2, 5, 7), 1e-9); +CheckAutoDiffJacobian(f, randn(2, 5, 7)+1i*randn(2, 5, 7), 1e-8); x = randn(3, 2, 1, 3); f = @(y) pagemtimes(x,y); -CheckAutoDiffJacobian(f, randn(2, 5, 7, 1), 1e-9); +CheckAutoDiffJacobian(f, randn(2, 5, 7, 1), 1e-8); y = randn(2, 4, 1); f = @(x) pagemtimes(x, y); -CheckAutoDiffJacobian(f, randn(3, 2, 1), 1e-9); +CheckAutoDiffJacobian(f, randn(3, 2, 1), 1e-8); x = randn(3, 2, 1); f = @(y) pagemtimes(x,y); -CheckAutoDiffJacobian(f, randn(2, 4, 1), 1e-9); +CheckAutoDiffJacobian(f, randn(2, 4, 1), 1e-8); x = randn(3, 2, 1); f = @(y) pagemtimes(x,y); -CheckAutoDiffJacobian(f, randn(2, 5, 5), 1e-9); +CheckAutoDiffJacobian(f, randn(2, 5, 5), 1e-8); x = randn(3, 2, 7); f = @(y) pagemtimes(x,y); -CheckAutoDiffJacobian(f, randn(2, 5, 7), 1e-9); +CheckAutoDiffJacobian(f, randn(2, 5, 7), 1e-8); y = randn(2, 5, 7, 2); f = @(x) pagemtimes(x, y); CheckAutoDiffJacobian(f, randn(3, 2, 7, 2), 1e-8); f = @(x) pagemtimes(x,x); -CheckAutoDiffJacobian(f, randn(3, 3, 5), 1e-9); +CheckAutoDiffJacobian(f, randn(3, 3, 5), 1e-8); f = @(x) norm(x); -CheckAutoDiffJacobian(f, rand(1, 3), 1e-9); -CheckAutoDiffJacobian(f, [-0.2818003 , 0.00971297, -0.00271337], 1e-9) -%CheckAutoDiffJacobian(f,rand(3,2),1e-9); uses svd, not coded yet +CheckAutoDiffJacobian(f, randn(3,1)+1i*randn(3,1), 1e-8); +f = @(x) norm(x,1); +CheckAutoDiffJacobian(f, [-0.2818003 , 0.00971297, -0.00271337], 1e-8) +%CheckAutoDiffJacobian(f,rand(3,2),1e-8); disp('norm-matrix OK') %uses svd, not yet OK + % testing repmat f = @(x) repmat(x, [3, 2]); -CheckAutoDiffJacobian(f, rand(2, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(2, 3), 1e-8); f = @(x) repmat(x, 1, 1, 10); -CheckAutoDiffJacobian(f, ones(3,3), 1e-9); +CheckAutoDiffJacobian(f, ones(3,3), 1e-8); f = @(x) x(:); -CheckAutoDiffJacobian(f, rand(2, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(2, 3), 1e-8); %testing compatible size multiplication (i.e. using broadcasting) f = @(x) x .* [3, 4, 2]; -CheckAutoDiffJacobian(f, rand(2, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(2, 3), 1e-8); f = @(x) [3, 4, 2] .* x; -CheckAutoDiffJacobian(f, rand(2, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(2, 3), 1e-8); f = @(x) x(1, :) .* x; -CheckAutoDiffJacobian(f, rand(2, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(2, 3), 1e-8); % various tests f = @(x) x'; -CheckAutoDiffJacobian(f, randn(2, 3), 1e-9); +CheckAutoDiffJacobian(f, randn(2,3), 1e-8); +CheckAutoDiffJacobian(f, rand(2,3)+1i*randn(2,3), 1e-8); f = @(x) abs(x); -CheckAutoDiffJacobian(f, randn(2, 3), 1e-9); +CheckAutoDiffJacobian(f, randn(2, 3), 1e-8); +CheckAutoDiffJacobian(f, rand(2,3)+1i*randn(2,3), 1e-8); f = @(x) sqrt(x); -CheckAutoDiffJacobian(f, rand(2, 3), 1e-8); +x = randn(2,3); x = x+0.1*sign(x); % avoid singularity at 0 +CheckAutoDiffJacobian(f, x, 1e-8); +z = randn(3)+1i*randn(3); z = z+0.1*sign(z); % avoid singularity at 0 +CheckAutoDiffJacobian(f, z, 1e-8); f = @(x) cos(x); -CheckAutoDiffJacobian(f, rand(2, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(2, 3), 1e-8); +CheckAutoDiffJacobian(f, randn(2,3)+1i*randn(2,3), 1e-7); f = @(x) sin(x); -CheckAutoDiffJacobian(f, rand(2, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(2, 3), 1e-8); +CheckAutoDiffJacobian(f, randn(2,3)+1i*randn(2,3), 1e-7); f = @(x) tan(x); -CheckAutoDiffJacobian(f, rand(2, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(2, 3), 1e-8); +CheckAutoDiffJacobian(f, randn(2,3)+1i*randn(2,3), 1e-7); +% AJR avoid magnified errors in steep gradients f = @(x) acos(x); -CheckAutoDiffJacobian(f, rand(2, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(2, 3)-0.3, 1e-8); +CheckAutoDiffJacobian(f, randn(2, 3)+1i*randn(2,3)-0.3, 1e-8); f = @(x) asin(x); -CheckAutoDiffJacobian(f, rand(2, 3), 1e-7); +CheckAutoDiffJacobian(f, rand(2, 3)-0.3, 1e-8); +CheckAutoDiffJacobian(f, randn(2, 3)+1i*randn(2,3)-0.3, 1e-8); f = @(x) atan(x); -CheckAutoDiffJacobian(f, rand(2, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(2, 3), 1e-8); +CheckAutoDiffJacobian(f, randn(2, 3)+1i*randn(2,3), 1e-7); f = @(x) exp(x); -CheckAutoDiffJacobian(f, rand(2, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(2, 3), 1e-8); +CheckAutoDiffJacobian(f, randn(2, 3)+1i*randn(2,3), 1e-8); f = @(x) log(x); -CheckAutoDiffJacobian(f, rand(2, 3)+0.1, 1e-9); +x = randn(2,3); x = x+0.1*sign(x); % avoid singularity at 0 +CheckAutoDiffJacobian(f, x, 1e-8); +z = randn(3)+1i*randn(3); z = z+0.1*sign(z); % avoid singularity +CheckAutoDiffJacobian(f, z, 1e-8); f = @(x) tanh(x); -CheckAutoDiffJacobian(f, rand(2, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(2, 3), 1e-8); +CheckAutoDiffJacobian(f, randn(2, 3)+1i*randn(2,3), 1e-7); f = @(x) conj(x); -CheckAutoDiffJacobian(f, rand(2, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(2, 3), 1e-8); +CheckAutoDiffJacobian(f, randn(2, 3)+1i*randn(2,3), 1e-7); t = rand(3, 3); f = @(x) cat(1, x, x*2, t); -CheckAutoDiffJacobian(f, rand(2, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(2, 3), 1e-8); f = @(x) repmat(x, [3, 4]); -CheckAutoDiffJacobian(f, rand(2, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(2, 3), 1e-8); f = @(x) diag(x); -CheckAutoDiffJacobian(f, rand(4, 1), 1e-9); -CheckAutoDiffJacobian(f, rand(4, 4), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 1), 1e-8); +CheckAutoDiffJacobian(f, rand(4, 4), 1e-8); f = @(x) diff(x, 1, 2); -CheckAutoDiffJacobian(f, rand(4, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 3), 1e-8); f = @(x) diff(x, 1, 1); -CheckAutoDiffJacobian(f, rand(4, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 3), 1e-8); +CheckAutoDiffJacobian(f, randn(4, 3)+1i*randn(4,3), 1e-8); f = @(x) diff(x, 1, 3); -CheckAutoDiffJacobian(f, rand(4, 3, 5, 2), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 3, 5, 2), 1e-8); f = @(x) x(:, end); -CheckAutoDiffJacobian(f, rand(4, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 3), 1e-8); f = @(x) x(end, :); -CheckAutoDiffJacobian(f, rand(4, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 3), 1e-8); f = @(x) x(2, :); -CheckAutoDiffJacobian(f, rand(4, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 3), 1e-8); f = @(x) max(x); -CheckAutoDiffJacobian(f, rand(4, 1), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 1), 1e-8); +CheckAutoDiffJacobian(f, rand(4, 1)+1i*rand(4,1), 1e-8); a = rand(4, 3); f = @(x) max(a, x); -CheckAutoDiffJacobian(f, rand(4, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 3), 1e-8); a = rand(4, 3); f = @(x) max(x, a); -CheckAutoDiffJacobian(f, rand(4, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 3), 1e-8); f = @(x) max(x); -CheckAutoDiffJacobian(f, rand(4, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 3), 1e-8); f = @(x) max(x, -x); -CheckAutoDiffJacobian(f, rand(4, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 3), 1e-8); f = @(x) min(x); -CheckAutoDiffJacobian(f, rand(4, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 3), 1e-8); f = @(x) min(x); -CheckAutoDiffJacobian(f, rand(4, 1), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 1), 1e-8); a = rand(4, 3); f = @(x) min(a, x); -CheckAutoDiffJacobian(f, rand(4, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 3), 1e-8); a = rand(4, 3); f = @(x) min(x, a); -CheckAutoDiffJacobian(f, rand(4, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 3), 1e-8); f = @(x) min(x, -x); -CheckAutoDiffJacobian(f, rand(4, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 3), 1e-8); f = @(x) x - x(1, 2); -CheckAutoDiffJacobian(f, rand(4, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 3), 1e-8); f = @(x) x - 3; -CheckAutoDiffJacobian(f, rand(4, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 3), 1e-8); f = @(x) 3 - x; -CheckAutoDiffJacobian(f, rand(4, 3), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 3), 1e-8); f = @(x) x^2; -CheckAutoDiffJacobian(f, 3, 1e-9); +CheckAutoDiffJacobian(f, randn(3), 1e-8); -f = @(x) x^3; -CheckAutoDiffJacobian(f, rand(3, 3), 1e-9); +f = @(x) x^5; +CheckAutoDiffJacobian(f, rand(3), 1e-8); f = @(x) x.^2; -CheckAutoDiffJacobian(f, rand(3, 2), 1e-9); - +CheckAutoDiffJacobian(f, rand(3, 2), 1e-8); f = @(x) x * 2.5; CheckAutoDiffJacobian(f, randn(2, 3), 1e-8); @@ -251,17 +284,19 @@ CheckAutoDiffJacobian(f, randn(2, 3), 1e-8); f = @(x) power(3,x); -CheckAutoDiffJacobian(f, randn(2, 3), 1e-9); +CheckAutoDiffJacobian(f, randn(2, 3), 1e-8); a = randn(2, 3); f = @(x) power(a,x); CheckAutoDiffJacobian(f, rand(2, 3), 1e-8); -a = rand(2, 3); +a = 1+rand(2, 3); % avoid singularity for small powers f = @(x) power(x,a); -CheckAutoDiffJacobian(f, rand(2, 3), 1e-7); +x = rand(2, 3); +CheckAutoDiffJacobian(f, x, 1e-8); -f = @(x) power(x,x*2); +% AJR, 7/2026: avoid singularity for small powers +f = @(x) power(x,1+x*2); CheckAutoDiffJacobian(f, rand(2, 3), 1e-7); % test matrix product @@ -291,75 +326,82 @@ f = @(x) inv(x); CheckAutoDiffJacobian(f, [[1,2,3];[3,1,2];[0,4,5]], 1e-6); +% AJR: generate random not-near-singular complex matrix +n = 5; [U,S,V]=svd(randn(n)+1i*randn(n)); +CheckAutoDiffJacobian(f, U*diag(cumsum(0.3+rand(1,n)))*V', 1e-6); + -f = @(x) x / x(2, 2); -CheckAutoDiffJacobian(f, rand(3, 2), 1e-6); +% AJR, 7/2026: avoid singular for small x22 +f = @(x) x / (1+x(2, 2)); +CheckAutoDiffJacobian(f, rand(3, 2), 1e-8); f = @(x) x / 3; -CheckAutoDiffJacobian(f, rand(3, 2), 1e-9); +CheckAutoDiffJacobian(f, rand(3, 2), 1e-8); -f = @(x) x ./ x(2, 2); -CheckAutoDiffJacobian(f, rand(3, 2), 1e-7); +% AJR, 7/2026: avoid singular for small x22 +f = @(x) x ./ (1+x(2, 2)); +CheckAutoDiffJacobian(f, rand(3, 2), 1e-8); f = @(x) x ./ 3; -CheckAutoDiffJacobian(f, rand(3, 2), 1e-9); +CheckAutoDiffJacobian(f, rand(3, 2), 1e-8); f = @(x) 3 ./ x; -CheckAutoDiffJacobian(f, randn(3, 2), 1e-4); +x = randn(3,2);% AJR: avoid singularity by +sign +CheckAutoDiffJacobian(f, x+sign(x), 1e-8); f = @(x) x .* abs(x); -CheckAutoDiffJacobian(f, randn(3, 3), 1e-9); +CheckAutoDiffJacobian(f, randn(3, 3), 1e-8); f = @(x) x .* x(2, 2); -CheckAutoDiffJacobian(f, rand(3, 2), 1e-9); +CheckAutoDiffJacobian(f, rand(3, 2), 1e-8); f = @(x) x + x(:, 1); -CheckAutoDiffJacobian(f, rand(3, 2), 1e-9); +CheckAutoDiffJacobian(f, rand(3, 2), 1e-8); f = @(x) reshape(x, 3, 2); -CheckAutoDiffJacobian(f, rand(3, 2), 1e-9); +CheckAutoDiffJacobian(f, rand(3, 2), 1e-8); f = @(x) sort(x); -CheckAutoDiffJacobian(f, rand(3, 2), 1e-9); +CheckAutoDiffJacobian(f, rand(3, 2), 1e-8); f = @(x) sort(x, 1); -CheckAutoDiffJacobian(f, rand(3, 2), 1e-9); +CheckAutoDiffJacobian(f, rand(3, 2), 1e-8); f = @(x) sort(x, 2); -CheckAutoDiffJacobian(f, rand(3, 2), 1e-9); +CheckAutoDiffJacobian(f, rand(3, 2), 1e-8); f = @(x) sort(x); -CheckAutoDiffJacobian(f, rand(3, 1), 1e-9); +CheckAutoDiffJacobian(f, rand(3, 1), 1e-8); f = @(x) x(3, :, :); -CheckAutoDiffJacobian(f, rand(3, 2, 4), 1e-9); +CheckAutoDiffJacobian(f, rand(3, 2, 4), 1e-8); f = @(x) x(:, 2, :); -CheckAutoDiffJacobian(f, rand(3, 2, 4), 1e-9); +CheckAutoDiffJacobian(f, rand(3, 2, 4), 1e-8); f = @(x) sum(x); -CheckAutoDiffJacobian(f, rand(3, 1), 1e-9) -CheckAutoDiffJacobian(f, rand(1, 3), 1e-9) +CheckAutoDiffJacobian(f, rand(3, 1), 1e-8) +CheckAutoDiffJacobian(f, rand(1, 3), 1e-8) f = @(x) sum(x, 2); -CheckAutoDiffJacobian(f, rand(3, 2), 1e-9); +CheckAutoDiffJacobian(f, rand(3, 2), 1e-8); f = @(x) cumsum(x, 2); -CheckAutoDiffJacobian(f, rand(3, 2, 4), 1e-9); +CheckAutoDiffJacobian(f, rand(3, 2, 4), 1e-8); f = @(x) cumsum(x, 2); -CheckAutoDiffJacobian(f, rand(3, 2), 1e-9); +CheckAutoDiffJacobian(f, rand(3, 2), 1e-8); f = @(x) cumsum(x); -CheckAutoDiffJacobian(f, rand(3, 2), 1e-9); +CheckAutoDiffJacobian(f, rand(3, 2), 1e-8); f = @(x) cumsum(x); -CheckAutoDiffJacobian(f, rand(3), 1e-9); +CheckAutoDiffJacobian(f, rand(3), 1e-8); f = @(x) mean(x, 2); -CheckAutoDiffJacobian(f, rand(3, 2, 4), 1e-9); +CheckAutoDiffJacobian(f, rand(3, 2, 4), 1e-8); f = @(x) mean(x, 1); -CheckAutoDiffJacobian(f, rand(3, 2, 4), 1e-9); +CheckAutoDiffJacobian(f, rand(3, 2, 4), 1e-8); f = @(x) mean(x) ; CheckAutoDiffJacobian(f, randn(1, 3), 1e-8); @@ -371,94 +413,73 @@ CheckAutoDiffJacobian(f, randn(3, 2, 4), 2e-9); t = randn(3, 2, 4); f = @(x) x .* t; -CheckAutoDiffJacobian(f, randn(3, 2, 4), 1e-9); +CheckAutoDiffJacobian(f, randn(3, 2, 4), 1e-8); f = @(x) t .* x; -CheckAutoDiffJacobian(f, randn(3, 2, 4), 1e-9); +CheckAutoDiffJacobian(f, randn(3, 2, 4), 1e-8); -f = @(x) eig(0.5*(x + x')); % we need to symetrize the matrix here so that the input to eig stays symetric when doing finite differences -t = randn(4, 4); -CheckAutoDiffJacobian(f, t, 1e-8); -f = @(x) selectKthOutput(2, 1, f, x); -t = randn(3, 3); -CheckAutoDiffJacobian(f, t+t', 1e-7); +f = @(x) stack2Fn('eig',x); +%disp('AJR: checking AD of symmetric arrays') +x=randn(4); x=x+x'; +CheckAutoDiffJacobian(f, x, 1e-8); +x=randn(3)+1i*randn(3); x=x+x'; +CheckAutoDiffJacobian(f, x, 1e-8); +%disp('AJR: checking AD of non-symmetric arrays') +CheckAutoDiffJacobian(f, randn(4), 1e-7); +CheckAutoDiffJacobian(f, randn(3)+1i*randn(3), 1e-7); f = @(x) x'; -CheckAutoDiffJacobian(f, randn(3, 3), 1e-9); +CheckAutoDiffJacobian(f, randn(3, 3), 1e-8); f = @(x) permute(x, [3, 1, 2]); -CheckAutoDiffJacobian(f, randn(3, 2, 4), 1e-9); +CheckAutoDiffJacobian(f, randn(3, 2, 4), 1e-8); f = @(x) x - 1; -CheckAutoDiffJacobian(f, randn(3, 2, 4), 1e-9); +CheckAutoDiffJacobian(f, randn(3, 2, 4), 1e-8); f = @(x) x + 1; -CheckAutoDiffJacobian(f, randn(3, 2, 4), 1e-9); +CheckAutoDiffJacobian(f, randn(3, 2, 4), 1e-8); f = @(x) [x, x * 2]; -CheckAutoDiffJacobian(f, randn(3, 2, 4), 1e-9); +CheckAutoDiffJacobian(f, randn(3, 2, 4), 1e-8); f = @(x) [x; x * 2]; -CheckAutoDiffJacobian(f, randn(3, 2, 4), 1e-9); +CheckAutoDiffJacobian(f, randn(3, 2, 4), 1e-8); f = @(x) det(x); -CheckAutoDiffJacobian(f, randn(2, 2), 1e-9); - -f = @(x) det(x); -CheckAutoDiffJacobian(f, randn(2, 2), 1e-9); +CheckAutoDiffJacobian(f, randn(2, 2), 1e-8); +%AJR, 7/2026: include complex value test f = @(x) det(x); -CheckAutoDiffJacobian(f, randn(3, 3), 1e-8) +CheckAutoDiffJacobian(f, 2*randn(3)+1i*randn(3), 1e-7 ); f = @(x) det(x); -CheckAutoDiffJacobian(f, randn(4, 4), 1e-9); - -f = @(x) sinh(x); -CheckAutoDiffJacobian(f, randn(2, 2), 1e-9); - -f = @(x) sinh(x); -CheckAutoDiffJacobian(f, randn(3, 3), 1e-9); +CheckAutoDiffJacobian(f, randn(4, 4), 1e-8); f = @(x) sinh(x); CheckAutoDiffJacobian(f, randn(4, 4), 1e-8); f = @(x) cosh(x); -CheckAutoDiffJacobian(f, randn(2, 2), 1e-9); - -f = @(x) cosh(x); -CheckAutoDiffJacobian(f, randn(3, 3), 1e-9); - -f = @(x) cosh(x); -CheckAutoDiffJacobian(f, randn(4, 4), 1e-9); - -f = @(x) asinh(x); -CheckAutoDiffJacobian(f, randn(2, 2), 1e-9); - -f = @(x) asinh(x); -CheckAutoDiffJacobian(f, randn(3, 3), 1e-9); +CheckAutoDiffJacobian(f, randn(4, 4), 1e-8); f = @(x) asinh(x); -CheckAutoDiffJacobian(f, randn(4, 4), 1e-9); - -f = @(x) acosh(x); -CheckAutoDiffJacobian(f, randn(2, 2), 1e-8); - -f = @(x) acosh(x); -CheckAutoDiffJacobian(f, rand(3, 3), 1e-8); +CheckAutoDiffJacobian(f, randn(4, 4), 1e-8); +z = 2*randn(3)+1i*randn(3); %AJR, 7/2026: include complex value test +CheckAutoDiffJacobian(f, z, 1e-8 ); +% AJR, 7/2026: the following used to lack testing the real cases of x>1 +% also adjust tolerance near the sqrts at x=+-1 f = @(x) acosh(x); -CheckAutoDiffJacobian(f, rand(4, 4), 1e-9); - -f = @(x) atanh(x); -CheckAutoDiffJacobian(f, rand(2, 2), 1e-9); +z = 2*randn(3)+1i*randn(3); +CheckAutoDiffJacobian(f, z, 1e-8/min(abs(z(:).^2-1)) ); +x = randn(3); +CheckAutoDiffJacobian(f, x, 1e-8/min(abs(x(:).^2-1)) ); +% AJR avoid magnified errors in steep gradients f = @(x) atanh(x); -CheckAutoDiffJacobian(f, rand(3, 3), 1e-9); - -f = @(x) atanh(x); -CheckAutoDiffJacobian(f, rand(4, 4), 1e-9); +CheckAutoDiffJacobian(f, rand(4, 4)-0.3, 1e-8); % some other tests @@ -475,5 +496,86 @@ CheckAutoDiffJacobian(f, x0, 1e-8); - +f = @(x) fft(x); +dim = randi([3 9],1,2); +X = randn(dim)+1i*randn(dim); +CheckAutoDiffJacobian(f, X, 1e-8); + +f = @(x) ifft(x); +dim = randi([3 9],1,2); +X = randn(dim)+1i*randn(dim); +CheckAutoDiffJacobian(f, X, 1e-8); + +% AJR: uses extra function to stack result U,S,V +% and adjust tolerance when close singular values +f = @(x) stack3Fn('svd',x); +X = randn(randi([3 9],1,2)); +CheckAutoDiffJacobian(f, X, 1e-8*(1+1/min(abs(diff(svd(X))))) ); +disp('AD.svd OK for real') +dim = randi([3 9],1,2); +X = randn(dim)+1i*randn(dim); +CheckAutoDiffJacobian(f, X, 1e-8*(1+1/min(abs(diff(svd(X))))) ); +disp('AD.svd OK for cmplx') + + + +disp('**** Now check on Analytic/holomorphic/...') +analyticFns={'transpose','sin','cos','tan','exp','log','sqrt','diff','sum','mean','cumsum','sinh','cosh','tanh','asin','acos','atan','asinh','acosh','atanh'} +for fn=analyticFns, nonAnalyticCheck(fn{1}); end; + +fn = 'inv', n=randi([3 5]); +if nonAnalytic(fn,randn(n)) + disp([fn '(real-matrix) not analytic']), end +if nonAnalytic(fn,randn(n)+1i*randn(n)) + disp([fn '(cmplx-matrix) not analytic']), end + +fn = @(X) stack2Fn('eig',X) +X=randn(randi([3 5])); X=X+X'; +if nonAnalytic(fn,X) + disp(['eig(sym-real-matrix) not analytic']), end +n=randi([2 4]); X=randn(n)+1i*randn(n); X=X+X'; +if nonAnalytic(fn,X) + disp(['eig(sym-cmplx-matrix) not analytic']), end +X=randn(randi([3 5])); +if nonAnalytic(fn,X) + disp(['eig(gen-real-matrix) not analytic']), end +n=randi([2 4]); X=randn(n)+1i*randn(n); +if nonAnalytic(fn,X) + disp(['eig(gen-cmplx-matrix) not analytic']), end + +disp('** Reaffirm non-analytic functions follow') +for fn={'imag','real','conj','abs','sign','norm','ctranspose'} + nonAnalyticCheck(fn{1}); +end; + + + +function ABC = stack3Fn(fn,X) + % AJR, 18 Jul 2026 + [m,n] = size(X); + [A,B,C] = feval(fn,X); + if strcmp(fn,'svd') % return only k columns, and + % ensure col.max-V is real-positive to be consistent with AD.svd() + k = 1:min(m,n); + [~,j] = max(abs(C(:,k))); + rot = conj( sign(C(n*(k-1)+j)) ); + ABC = [ A(:,k).*rot; B(k,k); C(:,k).*rot ]; + else ABC = [A;B;C]; + end;%if svd +end%function stack3Fn + +function VD = stack2Fn(fn,X) + % AJR, 18 Jul 2026 + [V,D] = feval(fn,X); + % always need to sort in a consistent order whether real + % or complex --- especially for cases when all real + % e-vals are perturbed to complex + d = diag(D); + [d,j] = sort(d,'ComparisonMethod','real'); + D = diag(d); V = V(:,j); + % ensure consistent normalised e-vecs, make max-abs element =1 + [~,iM]=max(abs(V)); + for j=1:size(V,2), V(:,j) = V(:,j)/V(iM(j),j); end + VD = [V D]; +end%function stack2Fn diff --git a/src/README.txt b/src/README.txt new file mode 100644 index 0000000..64098a6 --- /dev/null +++ b/src/README.txt @@ -0,0 +1,119 @@ +This version of auto-differentiation for Matlab, and maybe +Octave, is due to revisions by Tony Roberts, +http://orcid.org/0000-0001-8930-1552 + +* warning: some of these AutoDiff functions are not +analytic/holomorphic/... For those functions the +'perturbation variable' of the derivative is assumed to be +real. For those functions, to compute derivatives with +respect to a comples perturbation one must explicitly +include columns for pure-imaginary perturbations. For +example, z=AutoDiff(1+2i,[1 1i]), sign(z) cf sin(z) + +* extension of the base established by Martin de La Gorce, 2016, + +* with additions by Chris Noble, C2025, https://github.com/noblec04/MatlabGP + +* and revisions by Tony Roberts, 2026, to give this version +https://github.com/uoa1184615/EquationFreeGit/AutoDiff + +Sometimes info obtained by command such as "help AutoDiff.svd" + +2026-07-14 grep -e '^ *function ' AutoDiff.m +gives following list of AD functions: + function x = AutoDiff(values, derivatives) + function Jac = getderivs(x) + function val = getvalue(x) + function x = setdervis(x, derivatives) + function x = double(~) + function y = isinf(x) + function x = sinh(x) + function x = cosh(x) + function x = asinh(x) + function x = abs(x) + function x = acosh(x) + function x = atanh(x) + function x = sqrt(x) + function x = cos(x) + function x = sin(x) + function x = tan(x) + function x = acos(x) + function x = asin(x) + function y = ceil(x) + function y = floor(x) + function x = real(x) + function x = imag(x) + function x = atan(x) + function x = erf(x) + function x = erfc(x) + function x = gamma(x) + function x = exp(x) + function x = log(x) + function x = tanh(x) + function x = conj(x) + function y = fft(x,varargin) + function y = ifft(x,varargin) + function b = isreal(x) + function y = cat(dim, varargin) + function x = repmat(x, varargin) + function x = ctranspose(x) + function D = spdiags(B, d, m, n) + function D = diag(M) + function x = diff(x, n, dim) + function idx = end (x, k, n) + function z = eq(x, y) + function z = ne(x, y) + function z = sign(x) + function x = subsindex(x) + function z = ge(x, y) + function z = gt(x, y) + function z = le(x, y) + function z = lt(x, y) + function y = isnan(x) + function mylength = length(x) + function [m, id] = max(C, B) + function [m, id] = min(C, B) + function x = minus(x, y) + function x = mpower(x, n) + function x = inv(x) + function x = pinv(x,tol) + function [y,flag] = chol(x,shape) + function z = mldivide(x, y) + function z = mtimes(x, y) + function z = pagemtimes(x, y) + function z = mrdivide(x, y) + function x = norm(x, p) + function [Uad, Sad, Vad] = svd(A) + function n = numel(x) + function x = plus(x, y) + function x = power(x, y) + function x = rdivide(x, y) + function x = reshape(x, varargin) + function varargout = size(x, varargin) + function x = cumsum(x, varargin) + function varargout = sort(x, varargin) + function y = subsasgn(y, S, x) + function x = subsref(x, s) + function x = sum(x, dim) + function x = mean(x, dim) + function z = repmat_to_size(x, target_size) + function z = repmat_as(x, y) + function n = ndims(x) + function z = times(x, y) + function [V, D] = eig(C) + function x = transpose(x) + function x = permute(x, l) + function x = uminus(x) + function x = uplus(x) + function y = horzcat(varargin) + function y = det(x) + function y = vertcat(varargin) + function x = ones(varargin) + function x = zeros(varargin) + function x = nan(varargin) + function r = rank(x, varargin) + function C = pageMmult(A,At,B,Bt) + function M = spDiagFromVec(v) + function d = spdiag(a) + function D = transposeDiff(sizeM) + function M = subscriptDiff(idx, idy, sizeB) diff --git a/src/nonAnalytic.m b/src/nonAnalytic.m new file mode 100644 index 0000000..00b4369 --- /dev/null +++ b/src/nonAnalytic.m @@ -0,0 +1,24 @@ +function [flag] = nonAnalytic(fn,z) +% classifies whether a coded function fn is analytic or not +% at the complex point / vector / matrix z. Tests whether fn +% is analytic via auto-differentiation of the fn as provided +% by the package AutoDiff. Reports flag=0 if fn satisfies +% the Cauchy--Riemann eqns; =1 if fn fails u_x=v_y; =2 if +% fn fails u_y=-v_x; and =3 if fn fails both. +% AJR, 18 Jul 2026 +% Input: fn = string name or function handle of function +% z = complex/real-valued point/vector/matrix +% Output: flag = 0 if analytic ar z, otherwise 1,2,3 +[m,n] = size(z); +Z = AutoDiff( [real(z),imag(z)] ); +Z = Z(:,1:end/2)+1i*Z(:,end/2+1:end); +Fval = feval(fn,Z); +U = real(Fval); V=imag(Fval); +Du = reshape(full(getderivs(U)),[],m*n,2); +Dv = reshape(full(getderivs(V)),[],m*n,2); +uxMvy=Du(:,:,1)-Dv(:,:,2); +uyPvx=Du(:,:,2)+Dv(:,:,1); +uxMvyErr=norm(uxMvy,'fro'); +uyPvxErr=norm(uyPvx,'fro'); +flag = (uxMvyErr>1e-8)+2*(uyPvxErr>1e-8); +end%function nonAnalytic diff --git a/src/nonAnalyticCheck.m b/src/nonAnalyticCheck.m new file mode 100644 index 0000000..c19350c --- /dev/null +++ b/src/nonAnalyticCheck.m @@ -0,0 +1,21 @@ +function [flags] = nonAnalyticCheck(fn) +% classifies whether a coded function fn is analytic or not +% at random real/complex vector/matrix. Uses nonAnalytic() +% AJR, 17 Jul 2026 +% Input: fn = string name of function +% Output: flags = 0 if analytic ar z, otherwise 1,2,3 +m = randi([2 4]); +n = randi([3 5]); +flag = nonAnalytic(fn,randn(m,1)); +flags = flag; +if flag, disp([fn '(real-vector) not analytic ' num2str(flag)]), end +flag = nonAnalytic(fn,randn(m,1)+1i*randn(m,1)); +flags = [flags flag]; +if flag, disp([fn '(cmplx-vector) not analytic ' num2str(flag)]), end +flag = nonAnalytic(fn,randn(m,n)); +flags = [flags flag]; +if flag, disp([fn '(real-matrix) not analytic ' num2str(flag)]), end +flag = nonAnalytic(fn,randn(m,n)+1i*randn(m,n)); +flags = [flags flag]; +if flag, disp([fn '(cmplx-matrix) not analytic ' num2str(flag)]), end +end%function nonAnalyticCheck