Documentation for expm, logm and sqrtm

Author: mfasi

base/linalg/dense.jl

@@ -308,7 +308,7 @@ function logm(A::StridedMatrix)
         end
     end
     if np_real_eigs
-        warn("Matrix with nonpositive real eigenvalues, a nonprimary matrix logarithm will be returned.")
+        warn("Matrix with nonpositive real eigenvalues, a nonprincipal matrix logarithm will be returned.")
     end
 
     if isreal(A) && ~np_real_eigs

doc/stdlib/linalg.rst

@@ -638,17 +638,45 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
 .. function:: expm(A)
 
-   Compute the matrix exponential of ``A``.
+   Compute the matrix exponential of ``A``, defined by
+
+         .. math::
+
+            e^A = \sum_{n=0}^{\inf} \dfrac{A^n}{n!}.
+
+   For symmetric or Hermitian ``A``, an eigendecomposition (:func:`eigfact`) is used, otherwise the scaling and squaring algorithm (see [H05]_) is chosen.
+
+   .. [H05] Nicholas J. Higham, "The squaring and scaling method for the matrix
+      exponential revisited", SIAM Journal on Matrix Analysis and Applications,
+      26(4), 2005, 1179-1193.
+      `doi:/10.1137/090768539 <http://dx.doi.org/10.1137/090768539>`
 
 .. function:: logm(A)
 
-   Compute the matrix logarithm of ``A``.
+   If ``A`` has no negative real eigenvalue, compute the principal matrix logarithm of ``A``, i.e. the unique matrix :math:`X` such that :math:`e^X = A` and :math:`-\pi < Im(\lambda) < \pi` for all the eigenvalues :math:`\lambda` of :math:`X`. If ``A`` has nonpositive eigenvalues, a warning is printed and whenever possible a nonprimary matrix function is returned.
+
+   If `A` is symmetric or Hermitian, its eigendecomposition (:func:`eigfact`) is used, if `A` is triangular an improved version of the inverse scaling and squaring method is employed (see [AH12]_ and [AHR13]_). For general matrices, the complex Schur form (:func:`schur`) is computed and the triangular algorithm is used on the triangular factor.
+
+
+   .. [AH12] Awad H. Al-Mohy and Nicholas J. Higham, "Improved inverse  scaling
+      and squaring algorithms for the matrix logarithm", SIAM Journal on
+      Scientific Computing, 34(4), 2012, C153-C169.
+      `doi:/10.1137/110852553 <http://dx.doi.org/10.1137/110852553>`
+   .. [AHR13] Awad H. Al-Mohy, Nicholas J. Higham and Samuel D. Relton,
+      "Computing the Fréchet derivative of the matrix logarithm and estimating
+      the condition number", SIAM Journal on Scientific Computing, 35(4), 2013,
+      C394-C410.
+      `doi:/10.1137/120885991 <http://dx.doi.org/10.1137/120885991>`
 
 .. function:: sqrtm(A)
 
-   Compute the matrix square root of ``A``. If ``B = sqrtm(A)``, then ``B*B == A`` within roundoff error.
+   If ``A`` has no negative real eigenvalues, compute the principal matrix square root of ``A``, that is the unique matrix :math:`X` with eigenvalues having positive real part such that :math:`X^2 = A`. Otherwise, a nonprincipal square root is returned.
+
+   If `A` is symmetric or Hermitian, its eigendecomposition (:func:`eigfact`) is used to compute the square root. Otherwise, the square root is determined by means of the Björck-Hammarling method, which computes the complex Schur form (:func:`schur`) and then the complex square root of the triangular factor.
 
-   ``sqrtm`` uses a polyalgorithm, computing the matrix square root using Schur factorizations (:func:`schurfact`) unless it detects the matrix to be Hermitian or real symmetric, in which case it computes the matrix square root from an eigendecomposition (:func:`eigfact`). In the latter situation for positive definite matrices, the matrix square root has ``Real`` elements, otherwise it has ``Complex`` elements.
+  .. [BH83] Åke Björck and Sven Hammarling, "A Schur method for the square root
+     of a matrix", Linear Algebra and its Applications, 52-53, 1983, 127-140.
+     `doi:/10.1016/0024-3795(83)80010-X <http://dx.doi.org/10.1016/0024-3795(83)80010-X>`
 
 .. function:: lyap(A, C)
 

test/linalg3.jl

@@ -179,12 +179,12 @@ for elty in (Float64, Complex{Float64})
     A5  = convert(Matrix{elty}, [1 1 0 1; 0 1 1 0; 0 0 1 1; 1 0 0 1])
     @test_approx_eq expm(logm(A5)) A5
     s = readline(rd)
-    @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprimary matrix logarithm will be returned.")
+    @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprincipal matrix logarithm will be returned.")
 
     A6  = convert(Matrix{elty}, [-5 2 0 0 ; 1/2 -7 3 0; 0 1/3 -9 4; 0 0 1/4 -11])
     @test_approx_eq expm(logm(A6)) A6
     s = readline(rd)
-    @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprimary matrix logarithm will be returned.")
+    @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprincipal matrix logarithm will be returned.")
     redirect_stderr(OLD_STDERR)
 
     A7  = convert(Matrix{elty}, [1 0 0 1e-8; 0 1 0 0; 0 0 1 0; 0 0 0 1])