Why does MatrixFunction with Sinc return this error?

Question

The following:

mat = {{1., -3, -2}, {-1, 3, -2}, {0., 0., 0.}}
MatrixFunction[Sinc, mat]

returns:

MatrixFunction::drvnnum: Unable to compute the matrix function because 
the subsequent derivative of the function Sinc[#1] & at a numeric value is 
not numeric.

The documentation for this error gives examples where the function is not differentiable. But here, Sinc is, so why do I get this error?

Answer 1

Two things come into play:

1. Sinc'[0] yields Indeterminate, when it should be $0$.

2. This value needs to be evaluated in the computation of MatrixFunction[Sinc, mat].

In details: for the computation of MatrixFunction[] for an inexact matrix $\mathbf A$, a preliminary step is to always perform the Schur decomposition $\mathbf A=\mathbf Q\mathbf T\mathbf Q^\top$ (effectively, a similarity transformation into a triangular matrix). The reason for this is that the computation of a matrix function is easiest when a matrix is triangular.

Here, then is the required Schur decomposition:

MatrixForm /@ ({q, t} = SchurDecomposition[mat])

The problem arises in the computation of MatrixFunction[Sinc, t], since as mentioned previously, the evaluation of Sinc' at $0$ is needed, as evidenced by the following general formula:

MatrixFunction[f, {{0, a, b}, {0, c, d}, {0, 0, 0}}] // MatrixForm

A workaround in this case would be to use the explicit piecewise form of Sinc[]:

mySinc = Piecewise[{{1, # == 0}, {Sin[#]/#, True}}] &;
MatrixFunction[mySinc, mat]

You can check that the result of q.MatrixFunction[mySinc, t].Transpose[q] is the same.

However, this is not robust. Consider the following evaluation:

MatrixFunction[f, {{0, a, b}, {0, 0, d}, {0, 0, 0}}] // MatrixForm

where it can be seen that the second derivative is now involved. Neither Sinc'' nor mySinc'' will yield the right answer if evaluated on a matrix like this or similar to this, since one function is indeterminate and the other gives $0$ when the correct value should be $-\frac13$. For matrices of larger dimension where even higher derivatives might be involved, this is even more troublesome.

---

What to do when MatrixFunction[] fails?

I've always maintained that one of the nice things about Mathematica is that if all else fails, one can always go back to the mathematical definition. The situation is no different for MatrixFunction[]; as I have previously mentioned here and in other threads, the function seems to be based on evaluating a given function on a triangular matrix that is similar to the original matrix given to it, generated through either JordanDecomposition[] (in the exact case) or SchurDecomposition[] (in the inexact case). Since these seem to be failing here, we could use a different definition of the function of a matrix. Luckily, there are several of them.

In particular, one of the most useful definitions of a matrix function is the following Cauchy-like contour integral:

$$f(\mathbf A) = \frac{1}{2\pi i} \oint_\gamma f(z)\, (z \mathbf I- \mathbf A)^{-1}\,\mathrm dz$$

where $\gamma$ is a closed contour enclosing the eigenvalues of $\mathbf A$, and where $f(z)$ is analytic within. (I previously mentioned this here.)

As it turns out, this definition is useful for both exact and inexact computations. Let's deal with the exact case first, using an exact version of the matrix in the OP:

mat = {{1, -3, -2}, {-1, 3, -2}, {0, 0, 0}};

A slick way to use the contour integral definition is to recognize that the required computation is equivalent to the sum of the residues of the integrand at the eigenvalues of the original matrix, by virtue of the residue theorem and the Cauchy integral theorem.

Thus, we can compute the supposed result of MatrixFunction[Sinc, mat] in this way:

Sum[Map[Residue[#, {z, λ}] &, Sinc[z] Inverse[z IdentityMatrix[Length[mat]] - mat], {2}],
    {λ, Union[Eigenvalues[mat]]}]

(Note that Residue[] had to be mapped manually, as it does not automatically thread over lists.)

---

To demonstrate the application of the contour integral method to the inexact case, allow me to consider a different matrix and a different function with a removable singularity.

In particular, consider the matrix

rm = {{2/3, 0, -2/3, -1/3}, {0, 1, 0, 0}, {1/3, 0, 2/3, -2/3}, {2/3, 0, 1/3, 2/3}};

and suppose that you want to evaluate the function $f(\mathbf A)=(\mathbf A-\mathbf I)^{-1}\log\mathbf A$. The eigenvalues of this matrix are

eig = Eigenvalues[rm]
   {(-1)^(1/3), -(-1)^(2/3), 1, 1}

and thus, MatrixFunction[Log[#]/(# - 1) &, rm] is not expected to work in this case (and neither will Inverse[] + MatrixLog[]). For numerical evaluation, we use NIntegrate[] along with an appropriate choice of contour.

A particularly convenient contour can be obtained by treating the eigenvalues as points in the complex plane, and then slightly dilating the axis-aligned bounding rectangle that contains the eigenvalues. This can be done like so:

With[{ε = 1/20}, 
     corners = #1 + I #2 & @@@ Tuples[Transpose[(List @@ 
               BoundingRegion[ReIm[eig], "MinRectangle"]) + {-ε, ε}]][[{2, 1, 3, 4, 2}]]];

after which, one can then evaluate the contour integral:

NIntegrate[Log[z]/(z - 1) Inverse[z IdentityMatrix[4] - rm],
           {z, Sequence @@ corners} // Evaluate]/(2 π I) // Chop
   {{0.9379331214115499, 0, 0.271266454744782, 0.3333333333333591},
    {0, 1.0000000000001494, 0, 0},
    {-0.3333333333333591, 0, 0.9379331214115499, 0.271266454744782},
    {-0.271266454744782, 0, -0.3333333333333591, 0.9379331214115499}}

(Some NIntegrate::izero messages are thrown, but these are due to the zero entries and are harmless in this case.)

This compares favorably with the result of the residue-based method:

Sum[Map[Residue[#, {z, λ}] &, Log[#]/(# - 1) &[z] Inverse[z IdentityMatrix[4] - rm], {2}],
    {λ, Union[Eigenvalues[rm]]}] // FullSimplify

N[%]
   {{0.9379331214114057, 0., 0.2712664547447392, 0.3333333333333333},
    {0., 1., 0., 0.},
    {-0.3333333333333333, 0., 0.9379331214114057, 0.2712664547447392},
    {-0.2712664547447392, 0., -0.3333333333333333, 0.9379331214114057}}

---

For people who are interested in further details, I will recommend reading Functions of Matrices: Theory and Computation by Nick Higham.

Answer 2

Here's a workaround I've used for Sinc. It defines mySinc[k, x] to be the k-th derivative of Sinc[x], for k a nonnegative integer. (Update notice: Thanks to J.M. for pointing out a simple formula for the derivatives at zero.)

ClearAll[mySinc];
SetAttributes[mySinc, NHoldFirst];
mySinc[x_] := mySinc[0, x];
mySinc[0, x_?NumericQ] := Sinc[x];
mySinc[k_Integer?Positive, 0.] := N[(-1)^k Cos[k π/2]/(k + 1)];
mySinc[k_Integer?Positive, x_] /; x == 0 := 
  SetAccuracy[(-1)^k Cos[k π/2]/(k + 1),
   Piecewise[{{Log10[2 (k + 3)] + 2 Accuracy[x], EvenQ[k]}}, Log10[k + 2] + Accuracy[x]]];
mySinc[k_Integer?Positive, x_?NumericQ] /; x != 0 := 
  Evaluate@D[(# Cos[#] - Sin[#])/#^2, {#, k - 1}] &@x;
Derivative[n_Integer?Positive][mySinc][x_] := mySinc[n, x];
Derivative[0, n_Integer?Positive][mySinc][k_, x_] := mySinc[k + n, x];

Examples:

MatrixFunction[mySinc, mat] // MatrixForm

MatrixFunction[mySinc, Round@mat] // Simplify // MatrixForm

---

Caveat: The derivatives of Sinc are not numerically well-behaved near zero.

Labeled[
 TableForm[
  Table[SetPrecision[
    100 - Precision@Derivative[k][Sinc][SetPrecision[x, 100]],
    3], {x, 10^-Range[8]}, {k, 7}],
  TableHeadings -> {"10"^Row[{"-", #}] & /@ Range[8], Automatic}],
 {"x", "Derivative"}, {Left, Top}]

Precision loss of $\hbox{Sinc}^{(k)}(x)$ for $x$ near $0$.

The loss seems fairly predictable, and one could modify the definitions above to use sufficient precision internally to compute a result to a desired precision.

Answer 3

EDIT #2

The regularization method is more generally applicable than the Cauchy method.

Here's an example for which (according to J.M.) the Cauchy method fails but the regularization method works.

The matrix is that of the OP (integer form)

matx = {{1, -3, -2}, {-1, 3, -2}, {0, 0, 0}};

Eigenvalues[matx]

(* Out[90]= {4, 0, 0} *)

and the function has an essential singularity at 0. We obatin immediately

Limit[MatrixFunction[1/# Exp[-1/#] &, matx + \[Epsilon] IdentityMatrix[3]], \[Epsilon] -> 0] // FullSimplify

(* Out[50]= {{1/(16 E^(1/4)), -(3/(16 E^(1/4))), 1/(16 E^(1/4))}, {-(1/(16 E^(1/4))), 3/(16 E^(1/4)), -(1/(16 E^(1/4)))}, {0, 0, 0}} *)

Really simple, and fast.

EDIT

J.M.'s example treated with the regularization method

The Matrix:

rm = {{2/3, 0, -(2/3), -(1/3)}, {0, 1, 0, 0}, {1/3, 0, 2/3, -(2/3)}, {2/3, 0, 1/3, 2/3}}

The task is to calculate the matrix function

MatrixFunction[Log[#]/(# - 1) &, rm]

During evaluation of In[29]:= MatrixFunction::fnand: The function Log[#1]/(#1-1)& is not analytic or defined at 1. >>

The direct approach fails.

But regularization makes life easy:

MatrixFunction[Log[#]/(# - 1) &, 
  rm + \[Epsilon] IdentityMatrix[4]];

Limit[%, \[Epsilon] -> 0] // FullSimplify

(* Out[33]= {{1/9 (3 + Sqrt[3] \[Pi]), 0, 1/9 (-3 + Sqrt[3] \[Pi]), 1/
  3}, {0, 1, 0, 0}, {-(1/3), 0, 1/9 (3 + Sqrt[3] \[Pi]), 
  1/9 (-3 + Sqrt[3] \[Pi])}, {1/9 (3 - Sqrt[3] \[Pi]), 0, -(1/3), 
  1/9 (3 + Sqrt[3] \[Pi])}} *)

Summary

It may seem that the problem comes from the singularity of the Matrix mat, and is naturally related to the problem with Sin[x]/x/.x->0. But J.M. showed that problems arise if the matrix is defective.

There is no problem with MatrixFunction[Sinc, mat1] if Det[mat1] != 0.

We shall proceed here without exploring the root cause and show that the MatrixFunction of the OP can be computed exactly using two methods (a) power series, and (b) regularization, with the result:

MatrixFunction[Sinc, mat] == 

$$\begin{pmatrix} \frac{1}{16} (\sin (4)+12) & \frac{1}{16} (-3) (\sin (4)-4) & \frac{1}{16} (\sin (4)-4) \\ \frac{1}{16} (4-\sin (4)) & \frac{1}{16} (3 \sin (4)+4) & \frac{1}{16} (4-\sin (4)) \\ 0 & 0 & 1 \\ \end{pmatrix}$$

Derivation

Power series

We shall calculate the matrix function by a power series. This requires nothing but the powers of the matrix which always exist for a square matrix.

Here is the power series of Sinc[]:

Sum[(-1)^k x^(2 k)/(2 k + 1)!, {k, 0, ∞}]

(* Out[10]= Sin[x]/x *)

The matrix in question is

mat = {{1., -3, -2}, {-1, 3, -2}, {0.`, 0.`, 0.`}} 

(* Out[9]= {{1., -3, -2}, {-1, 3, -2}, {0., 0., 0.}} *)

In order to calculate the powers symbolically we take numerical constants and define

mat1 = Floor[mat]

(* Out[11]= {{1, -3, -2}, {-1, 3, -2}, {0, 0, 0}} *)

From the first few even powers we can easily guess the general form:

mp[k_] := {{4^(2 k - 1), -3 4^(2 k - 1), 4^(2 k - 1)}, {-4^(2 k - 1), 
    3 4^(2 k - 1), -4^(2 k - 1)}, {0, 0, 0}};

Notice that, due to the singularity of mat1, the zeroth power does not exist.

MatrixPower[mat1, 0]

During evaluation of In[15]:= MatrixPower::sing: Matrix {{1,-3,-2},{-1,3,-2},{0,0,0}} is singular. >>

(* Out[15]= MatrixPower[{{1, -3, -2}, {-1, 3, -2}, {0, 0, 0}}, 0] *)

Hence we take care of the zeroth power by simply adding the unit matrix. We find finally

f = DiagonalMatrix[Array[1 &, 3]] + 
  Sum[(-1)^k mp[k]/(2 k + 1)!, {k, 1, ∞}]

(* Out[46]= {{1 + 1/16 (-4 + Sin[4]), -(3/16) (-4 + Sin[4]), 
  1/16 (-4 + Sin[4])}, {1/16 (4 - Sin[4]), 1 + 3/16 (-4 + Sin[4]), 
  1/16 (4 - Sin[4])}, {0, 0, 1}} *)

Regularization

If we make the matrix into a regular one by adding a small matrix

mat2 = mat1 + ε DiagonalMatrix[Array[1 &, 3]];

Det[mat2]

(* Out[48]= ε (4 ε + ε^2) *)

we can apply the MatrixFunction without any problem:

f2 = MatrixFunction[Sinc, mat2];

and then take the limit

f20 = Limit[f2, ε -> 0];

This gives the same result as before

f == f20

(* Out[52]= True *)

Numerical matrix, partial series

Defining the partial series as

fn[nn_] := 
 DiagonalMatrix[Array[1 &, 3]] + 
  Sum[(-1)^k MatrixPower[mat, 2 k]/(2 k + 1)!, {k, 1, nn}]

we can calculate the numerical result up to the nn-th term.

For nn = 8 the result is alrady we find

fn[8] // N

(* Out[64]= {{0.7027, 0.8919, -0.2973}, {0.2973, 0.1081, 0.2973}, {0., 0., 1.}} *)

in fair agreement with the exact values:

f // N

(* Out[65]= {{0.7027, 0.8919, -0.2973}, {0.2973, 0.1081, 0.2973}, {0., 0., 1.}} *)

Numerical matrix, regularization

Regularizing the numerical matrix of the OP gives

mat4 = mat + ε DiagonalMatrix[Array[1 &, 3]];

Application of MatrixFunction with given small ε

f4 = MatrixFunction[Sinc, mat4] /. ε -> 10^(-8)

During evaluation of In[96]:= JordanDecomposition::jdimp: Unable to find the Jordan decomposition of the matrix with the given precision. Try higher precision or SchurDecomposition instead. >>

(* Out[96]= {{0.7027, 0.8919, -0.2973}, {0.2973, 0.1081, 0.2973}, {0., 0., 1.}} *)

results in the known numerical result despite of the error message.

Examples of singular matrices

An extreme case of a matrix is

mat3 = {{0, 0}, {0, 0}};

Nevertheless, there's no problem with Sinc:

MatrixFunction[Sinc, mat3]

(* Out[73]= {{1, 0}, {0, 1}} *)

Matrices of the type {{a,b},{0,0}} with a,b [Element]{0,1}

mm = {#, {0, 0}} & /@ Union[Tuples[{1, 0, 0, 0}, 2]]

(* Out[15]= {{{0, 0}, {0, 0}}, {{0, 1}, {0, 0}}, {{1, 0}, {0, 0}}, {{1, 1}, {0, 0}}} *)

Table[{k, MatrixFunction[Sinc, mm[[k]]]}, {k, 1, 4}]

During evaluation of In[18]:= MatrixFunction::fnanc: The function Sinc[#1]& is not analytic at 0. >>

(* Out[18]= {
{1, {{1, 0}, {0, 1}}}, 
{2, MatrixFunction[inc, {{0, 1}, {0, 0}}]}, 
{3, {{Sinc[1], 0}, {0, 1}}}, 
{4, {{Sinc[1], -1 + Sinc[1]}, {0, 1}}}
} *)

Notice that only for mm[[2]] Mathematica complains the presumed non-analyticity of Sinc at 0 and Returns the input unevaluated.

Again the regularization method removes the problem:

meps = \[Epsilon] IdentityMatrix[2];

Table[{k, MatrixFunction[Sinc, meps + mm[[k]]]}, {k, 1, 4}];

Limit[%, \[Epsilon] -> 0]

(* Out[21]= {{1, {{1, 0}, {0, 1}}}, {2, {{1, 0}, {0, 1}}}, {3, {{Sin[1], 0}, {0, 1}}}, {4, {{Sin[1], -1 + Sin[1]}, {0, 1}}}} *)