# Why does MatrixFunction with Sinc return this error?

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?

• Perhaps Sinc'[0.] and Sinc'[x] give a clue? (The derivative of Sinc unfortunately reverts to the derivative of Sin[x]/x.) – Michael E2 Oct 11 '16 at 2:33
• @MichaelE2 I see. But why does MatrixFunction requires the derivatives, and not "just" the sum of the ponderated matrix powers? – anderstood Oct 11 '16 at 2:47
• – Szabolcs Oct 11 '16 at 11:51
• @anderstood Good question. Please see my contribution. – Dr. Wolfgang Hintze Oct 11 '16 at 13:05
• I wish to add a comment but I do not have sufficient reputation to do so. The MATLAB function wkm can be used to compute sinc(A) as follows: A = [1, -3, -2; -1, 3, -2; 0, 0, 0]; [C, S] = wkm(-A^2); The matrices C and S returned by wkm are the matrix functions cos(A) and sinc(A). – npras Jun 21 '18 at 20:45

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.

• Would rationalizing be a workaround? – anderstood Oct 11 '16 at 3:54
• In that case, the Jordan decomposition is used instead of Schur, but derivatives will still be needed. – J. M.'s technical difficulties Oct 11 '16 at 3:59
• Well, I wish I could upvote again! – Michael E2 Oct 13 '16 at 1:01
• The note of appreciation suffices, @Michael. :) – J. M.'s technical difficulties Oct 13 '16 at 1:03
• An advantage of the contour integral definition that I forgot to mention is that both versions are easily modified to implement the "action form" of a matrix function, $f(\mathbf A)\mathbf x$: one merely replaces the Inverse[] with an appropriate call to LinearSolve[] on the given vector. – J. M.'s technical difficulties Oct 13 '16 at 1:12

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.

• Here's a slightly simpler expression: mySinc[k_Integer?Positive, 0.] := N[(-1)^k Cos[k π/2]/(k + 1)] – J. M.'s technical difficulties Oct 11 '16 at 14:44
• @J.M. Thanks! I did think of something like that, but I didn't think of using in mySinc[.., 0.]. Strange, them blind spots. – Michael E2 Oct 11 '16 at 14:53
• You might also be interested in this explicit formula: mySinc[k_Integer?Positive, x_?NumericQ] /; x != 0 := Sum[(-1)^j x^(-j - 1) j! Binomial[k, j] Sin[x + (k - j) π/2], {j, 0, k}]. – J. M.'s technical difficulties Oct 12 '16 at 11:16

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}}}} *)

• "from the singularity of the Matrix" - it's not necessarily due to singularity. For example, MatrixFunction[Sinc, {{2, 2, 3}, {0, 0, 1}, {0, 0, 1}}] evaluates without a hitch. It's when the matrix is defective that we get into trouble with Sinc's derivatives. – J. M.'s technical difficulties Oct 11 '16 at 11:28
• @J.M. Good counterexample to my first statement. Do you also have one for the case Det[mat]!=0 an still difficulties with Sinc? I think it is always possible to calculate the matrix function by a power series, if applicable. That's what I did here to find the result for the specific mat. – Dr. Wolfgang Hintze Oct 11 '16 at 11:40
• If there are no zero eigenvalues (thus, nonsingular), then there isn't any trouble with Sinc[] and its derivatives. Yes, power series are one way to evaluate matrix functions when the Schur-Parlett algorithm is not usable, but there are also other ways. – J. M.'s technical difficulties Oct 11 '16 at 11:44
• BTW: IdentityMatrix[] is built-in, so you can do IdentityMatrix[3] instead of DiagonalMatrix[Array[1 &, 3]]. For your last example, I would consider {{0, 1}, {0, 0}} to be the simplest example of a matrix whose Sinc[]` can't be evaluated. – J. M.'s technical difficulties Oct 13 '16 at 0:55
• @J.M. Thanks. (1) IdentityMatrix[] I once knew it ... -) (2) I have added some more examples, including yours, and also showed that regularization removes the analyticity problem with it. – Dr. Wolfgang Hintze Oct 13 '16 at 9:45