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