A recently developed numerical method for the calculation of derivatives of functions of general complex matrices, which can also be combined with implicit matrix function approximations such as Krylov-Ritz type algorithms, is presented. An important use case for the method in the context of lattice gauge theory is the overlap Dirac operator at finite quark chemical potential. Derivatives of the lattice Dirac operator are necessary for the computation of conserved lattice currents or the fermionic force in Hybrid Monte-Carlo and Langevin simulations. To calculate the overlap Dirac operator at finite chemical potential the product of the sign function of a non-Hermitian matrix with a vector has to be computed. For non-Hermitian matrices it is not possible to efficiently approximate the sign function with polynomials or rational functions. Implicit approximation algorithms, like Krylov-Ritz methods, that depend on the source vector have to be used instead. Our method can also provide derivatives of such implicit approximations. We show how a generalised deflation prescription can be used to improve the performance of the method, if some eigenvalues and eigenvectors of the matrix being differentiated are known. To show that the method is efficient and well suited for practical calculations we provide test results for the two-sided Lanczos approximation of the finite-density overlap Dirac operator on SU(3) gauge field configurations on lattices with sizes up to 14 × 14 3 .