Ygtjki

I have a matrix $A$ for which I want to compute the quantity $Tlambda_j = Pi_lambda_ine lambda_j fracA - lambda_i Ilambda_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]

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.

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_]]

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.