# Skipping indices in a product

I have a matrix $$A$$ for which I want to compute the quantity $$T\lambda_j = \Pi_{\lambda_i\ne \lambda_j} \frac{A - \lambda_i I}{\lambda_j-\lambda_i}$$, where $$\lambda_i$$ ($$\lambda_j$$) denote the eigenvalues of $$A$$. How can this be implemented in Mathematica? Just gave a try here:

A = {{1, 0, 0, 1},{0, 1, 2, 0},{1, 1, 0, 2},{0, 0, 0, 1}};
Eigenvalues[A]


{2, -1, 1, 1}

Tj = Product[(A - Eigenvalues[A][[i]] IdentityMatrix[4])/(
Eigenvalues[A][[j]] - Eigenvalues[A][[i]]), {i, 1, 4}]

• in which part exactly you want to exclude it in Tj !?? – Alrubaie Mar 29 '19 at 18:02
• do you want it to be skipped put not Zero right !? – Alrubaie Mar 29 '19 at 18:07
• @Alrubaie, there was a typo in my post. Just edited it. I want the denominator to be non-zero and hence avoid the case for which $i=j$. – H. Kenan Mar 29 '19 at 18:08
• @Alrubaie, my $i$ and $j$ are not the indices in my question. They are the eigenvalues. I should have used something like $\lambda_i$ and $\lambda_j$. – H. Kenan Mar 29 '19 at 18:12
• That product is presumably a matrix multiplication? – J. M.'s discontentment Mar 29 '19 at 18:34

Here is my pedestrian implementation of your formula:

a = {{1, 0, 0, 1}, {0, 1, 2, 0}, {1, 1, 0, 2}, {0, 0, 0, 1}};

ClearAll[t]
t[amat_, j_] := Module[
{evals, usable},
evals = Eigenvalues[amat];
usable = DeleteDuplicates@Cases[evals, Except@evals[[j]] ];
Dot @@
Table[
(amat - i IdentityMatrix[Length[amat]])/(evals[[j]] - i),
{i, usable}
]
]

t[a, 4]


You do not provide an example of desired output, so I will let you check whether this is what you expect.

• Thanks, @MarcoB. It leads precisely to the expected result. However, it looks too complicated. Nevertheless, it is fine as it works. – H. Kenan Mar 29 '19 at 18:46

Something like this?

Clear[A, evals, T]
A = {{1, 0, 0, 1}, {0, 1, 2, 0}, {1, 1, 0, 2}, {0, 0, 0, 1}};
T[A_?MatrixQ, j_Integer] := With[
{evals = Eigenvalues[A], id = IdentityMatrix@Length@A},
Dot @@ Table[
If[evals[[j]] - evals[[i]] == 0, id, (A - evals[[i]] id)/(evals[[j]] - evals[[i]])],
{i, Length@A}
]
]

MatrixForm /@ Array[T[A, #] &, 4]


This

A = {{1, 0, 0, 1},{0, 1, 2, 0},{1, 1, 0, 2},{0, 0, 0, 1}};
e=Eigenvalues[A];
Map[(A-e[[#[[1]]]]*IdentityMatrix[4])/(e[[#[[2]]]]-e[[#[[1]]]])&,
DeleteCases[Tuples[Range[4],2],{i_,i_}]]


generates your twelve matricies with i not equal to j.

Put Dot@@ in front of that Map to form the dot product of the 12 matricies.

That works by forming every possible distinct i,j pair and then using those in the Map

If it might be easier to read you can also write it this way

Map[(ei=e[[#[[1]]]];ej=e[[#[[2]]]];
(A-ei*IdentityMatrix[4])/(ej-ei))&,
DeleteCases[Tuples[Range[4],2],{i_,i_}]]

• Should e[[#[[2]]]]-e[[[[1]]]] be e[[#[[2]]]]-e[[#[[1]]]]? – NonDairyNeutrino Mar 29 '19 at 19:12
• @ThatGravityGuy Yes! Good catch. Thank you! Corrected. – Bill Mar 29 '19 at 19:14

Another way:

ClearAll[t];
t[j_Integer, A_?SquareMatrixQ] := t[j, A, Eigenvalues@A];   (* add the eigenvalues *)
t[j_Integer, A_?SquareMatrixQ, evs_?VectorQ] /; Length@A == Length@evs := (* arg checks *)=
Fold[
#1.(A - #2 IdentityMatrix[Length@A])/(evs[[j]] - #2) &,
IdentityMatrix[Length@A],
Pick[evs, Unitize[evs - evs[[j]]], 1]                    (* Pick nonzero differences *)
];


Performance tuning: One can use DeleteCases[evs, e_ /; e == evs[[j]]] to pick the eigenvalues that give a nonzero difference. It makes no consistent difference to the timing on a 101 x 101 machine real matrix. One can save a little time by computing the identity matrix once and using With[] to inject it in the two places it occurs. One can also save time using dot = (dot = Dot; #2) & instead of Dot to skip the multiplication by the identity matrix (the first step of Fold[]). The differences evs - evs[[j]] appear twice, so they can be replaced by a single computation like the identity matrix. It can make up to a 10% improvement.