dfunm

 Fréchet derivative of matrix functions.

 function [D, fX] = dfunm(funm, X, H)

 Computes the directional derivative (the Fréchet derivative) of a matrix
 function (such as @logm, @expm, ...) at X along H (square matrices),
 according to a very nice trick which appears in this paper:
 
 "Computing the Fréchet derivative of the matrix exponential, with an
 application to condition number estimation",
 Awad H. Al-Mohy and Nicholas J. Higham, 2009.
 http://eprints.ma.man.ac.uk/1218/01/covered/MIMS_ep2008_26.pdf

 Thus, D = lim_(t -> 0) (funm(X + tH) - funm(X)) / t.

 The second output is fX = funm(X), which comes out as a by-product. It
 may be less accurate than calling funm(X) directly.

 Note: under mild conditions, the adjoint of dfunm(X, .) is dfunm(X', .),
 which is a fact often useful to derive gradients of matrix functions
 involving funm(X).
 (This is wrt the inner product inner = @(A, B) real(trace(A'*B))).

 This code is simple, but may not be the most efficient. In particular, it
 requires computing the matrix function on matrices which are four times
 as big, and which may have lost important structure (such as symmetry).
 
 See also: dlogm dexpm dsqrtm

0001 function [D, fX] = dfunm(funm, X, H)
0002 % Fréchet derivative of matrix functions.
0003 %
0004 % function [D, fX] = dfunm(funm, X, H)
0005 %
0006 % Computes the directional derivative (the Fréchet derivative) of a matrix
0007 % function (such as @logm, @expm, ...) at X along H (square matrices),
0008 % according to a very nice trick which appears in this paper:
0009 %
0010 % "Computing the Fréchet derivative of the matrix exponential, with an
0011 % application to condition number estimation",
0012 % Awad H. Al-Mohy and Nicholas J. Higham, 2009.
0013 % http://eprints.ma.man.ac.uk/1218/01/covered/MIMS_ep2008_26.pdf
0014 %
0015 % Thus, D = lim_(t -> 0) (funm(X + tH) - funm(X)) / t.
0016 %
0017 % The second output is fX = funm(X), which comes out as a by-product. It
0018 % may be less accurate than calling funm(X) directly.
0019 %
0020 % Note: under mild conditions, the adjoint of dfunm(X, .) is dfunm(X', .),
0021 % which is a fact often useful to derive gradients of matrix functions
0022 % involving funm(X).
0023 % (This is wrt the inner product inner = @(A, B) real(trace(A'*B))).
0024 %
0025 % This code is simple, but may not be the most efficient. In particular, it
0026 % requires computing the matrix function on matrices which are four times
0027 % as big, and which may have lost important structure (such as symmetry).
0028 %
0029 % See also: dlogm dexpm dsqrtm
0030 
0031 % This file is part of Manopt: www.manopt.org.
0032 % Original author: Nicolas Boumal, July 3, 2015.
0033 % Contributors:
0034 % Change log:
0035 %
0036 %   June 14, 2019 (NB): now also outputs funm(X) as a by-product.
0037     
0038     n = size(X, 1);
0039     
0040     assert(length(size(X)) == 2,     'X and H must be square matrices.');
0041     assert(length(size(H)) == 2,     'X and H must be square matrices.');
0042     assert(size(X, 1) == size(X, 2), 'X and H must be square matrices.');
0043     assert(all(size(X) == size(H)),  'X and H must have the same size.');
0044     
0045     Z = zeros(n);
0046     A = funm([X, H ; Z, X]);
0047     D = A(1:n, (n+1):end);
0048     fX = A(1:n, 1:n);
0049     
0050 end