1------------------------------------------------------------
2-- TUTORIAL 2
3------------------------------------------------------------
4-- COMMUTING FAMILIES OF ENDOMORPHISMS
5------------------------------------------------------------
6-- ring for minimal polynomials
7QQz ::= QQ[z];
8------------------------------------------------------------
9-- EXAMPLE  our first commuting family of endomorphisms
10
11M := IdentityMat(QQ, 9);  M;
12L := GetRows(M);  L;
13
14Z3 := ZeroMat(QQ, 3, 3);
15
16A := [BlockMat([[Z3, MakeMatByCols(3,3,L[i])],
17	        [Z3,            Z3          ]])
18      | i in 1..9];
19
20A[1];
21A[2]; -- and so on..
22A[4]*A[7];  --> 0
23-- All products of two matrices gives 0-matrix
24--    therefore they trivially commute ;-)
25
26--------------------------------------------------
27-- DEFINITION
28-- Given a set of commuting endomorphisms S
29--   the COMMUTING FAMILY  F  GENERATED BY  S
30--     is the algebra F = K[S].
31--   It may be presented as K[z_1,..,z_n]/J
32--------------------------------------------------
33
34-- The QQ-basis of the family F is Id(6)=A[i]^0, A[1], ...A[9]
35-- therefore the dimension is 10
36
37-- F = K[A] ~ K[z_1,..,z_9]/J
38use QQ[z[1..9]];
39BM := BMForMat(A); -- generalization of MinPoly to a set of matrices:
40-- Groebner Basis (of J, the ideal of relations)
41-- Quotient Basis (basis as K-vector space of K[z_1,..,z_9]/J)
42indent(BM);
43
44len(BM.QB); --> 10
45Rel := BM.GB;
46-- evaluation of these polys in the matrices A gives 0-matrix
47
48------------------------------------------------------------
49-- RECALL from Tutorial1:
50define TutEvalPolyInMat(f, M)
51  ResM := ZeroMat(RingOf(M), NumRows(M), NumCols(M));
52  foreach summand in monomials(f) do
53    ResM := ResM + LC(summand) * M^deg(LPP(summand));
54  endforeach;
55  return ResM;
56enddefine; -- TutEvalPolyInMat
57
58------------------------------------------------------------
59-- NEW FUNCTION: eval multivariate poly in list of matrices
60--   (generalization of TutEvalPolyInMat)
61define TutEvalPolyInMatList(f, M)
62  ResM := ZeroMat(RingOf(M[1]), NumRows(M[1]), NumCols(M[1]));
63  foreach summand in monomials(f) do
64    d := exponents(LPP(summand));
65    ResM := ResM
66          + LC(summand) * product([ M[i]^d[i] | i in 1..len(M)]);
67  endforeach;
68  return ResM;
69enddefine; -- TutEvalPolyInMatList
70------------------------------------------------------------
71
72-- verify: some relations on A, the list of our 9 matrices
73IsZero( TutEvalPolyInMatList(Rel[1], A) );  --> true
74IsZero( TutEvalPolyInMatList(Rel[13], A) ); --> true
75
76------------------------------------------------------------
77-- THE KERNEL OF AN IDEAL IN A COMMUTING FAMILY
78
79--------------------------------------------------
80-- DEFINITION   I = <psi_1, ... psi_n> ideal in F.
81--   Ker(I) := intersection_i (Ker(psi_i))
82--------------------------------------------------
83
84M1 := BlockMat([[Z3, MakeMatByCols(3,3, [1,2,3, 4,5,6, 7,8,9])],
85		[Z3, Z3]]);
86M2 := BlockMat([[Z3, MakeMatByCols(3,3, [9,8,7, 6,5,4, 3,2,1])],
87		[Z3, Z3]]);
88M3 := BlockMat([[Z3, MakeMatByCols(3,3, [1,3,1, 1,3,1, 1,3,1])],
89		[Z3, Z3]]);
90C := BlockMat([[M1], [M2], [M3]]);   C;
91
92LinKerBasis(M1); -- basis of Ker(phi1)
93LinKerBasis(M2); -- basis of Ker(phi2)
94LinKerBasis(M3); -- basis of Ker(phi3)
95LinKerBasis(C);  -- basis of intersection of the 3 kernels
96
97-- EXERCISE: try with other endomorphisms (matrices)
98
99--------------------------------------------------
100-- ALGORITHM   1-DIMENSIONAL JOINT EIGENSPACES
101--------------------------------------------------
102
103-- input
104A1:=matrix(QQ,
105 [[0, 1, 0, 0],
106  [0, 0, 0, 1],
107  [-6, 5, 2, -1],
108  [4, -8, 0, 5]]);
109A2:= matrix(QQ,
110 [[0, 0, 1, 0],
111  [-6, 5, 2, -1],
112  [3, -4, 2, 1],
113  [-16, 12, 4, -2]]);
114A1*A2 = A2*A1;
115
116-- (1)
117use QQz;
118p1 := eigenfactors(A1, z);   p1;  -- 2 linear factors
119B1_1 := LinKer(TutEvalPolyInMat(p1[1], A1));
120B1_2 := LinKer(TutEvalPolyInMat(p1[2], A1));
121
122B1_1; B1_2;
123
124-- (2)
125S := [B1_1]; -- bases of length 1
126L := [B1_2]; -- bases of length >1
127
128-- loop on Ai, loop on bases in S
129
130-- (5)
131MatB := B1_2;
132Restr := A2 * MatB;  Restr;  -- restriction of phi2 to <B1>
133B1_2;
134MatRestr := LinSolve(MatB, Restr);
135MatRestr;
136
137-- (6)
138p21 := eigenfactors(MatRestr, z);  p21;
139LinKerBasis(TutEvalPolyInMat(p21[1], MatRestr));
140LK := LinKer(TutEvalPolyInMat(p21[1], MatRestr)); -- directly as matrix
141LK;
142JointEig2 := MatB * LK;
143JointEig2;
144
145-- (7)
146append(ref S, JointEig2);
147indent(S);
148
149-- verification:
150GetCol(S[1],1);
151GetCol(A1*S[1], 1);
152GetCol(A2*S[1], 1);
153LinSolve(S[1], A1*S[1]); -- eigenvalue of phi1
154LinSolve(S[1], A2*S[1]); -- eigenvalue of phi2
155
156GetCol(S[2],1);
157GetCol(A1*S[2], 1);
158GetCol(A2*S[2], 1);
159LinSolve(S[2], A1*S[2]); -- eigenvalue of phi1
160LinSolve(S[2], A2*S[2]); -- eigenvalue of phi2
161
162------------------------------------------------------------
163-- EXAMPLE
164
165A1 := mat([[0,2,0,0],
166	   [1,0,0,0],
167	   [0,0,0,2],
168	   [0,0,1,0]]);
169
170A2 := mat([[0,0,8,0],
171	   [0,0,0,8],
172	   [1,0,0,0],
173	   [0,1,0,0]]);
174A1*A2 = A2*A1;
175-- F commuting family generated by A1, A2
176
177use QQz;
178eigenfactors(A1, z);
179eigenfactors(A2, z);
180A3 := A1 + A2;  -- psi = phi1 + phi2
181eigenfactors(A3, z); --> two eigenfactors
182--> so phi1 and phi2 are not splitting endomorphisms
183
184p := eigenfactors(A3, z);
185NumCols(LinKer(TutEvalPolyInMat(p[1], A3))); --> dim 2
186NumCols(LinKer(TutEvalPolyInMat(p[2], A3))); --> dim 2
187--> 2+2 = dim(QQ^4)  --> so BigKer = LinKer
188
189-- EXERCISE check with the algorithm described above that
190--  there are no 1-dimensional joint eigenspaces of F
191------------------------------------------------------------
192-- ALGORITHM CYCLICITY TEST
193A1 := mat([[0,1,0,0],
194	   [0,0,0,1],
195	   [-1,1,1,0],
196	   [9,-12,-1,6]]);
197
198A2 := mat([[0,0,1,0],
199	   [-1,1,1,0],
200	   [-2,0,3,0],
201	   [-1,0,1,1]]);
202
203A1*A2 = A2*A1;
204-- F commuting family generated by A1, A2
205
206-- (1)
207RelRing ::= QQ[x,y];
208Rel := BMForMat(RelRing, [A1,A2]);
209indent(Rel);
210QB := Rel.QB;
211
212-- (3)
213AuxRing ::= QQ[w[1..4]];  -- auxiliary indeterminates
214W := ColMat(indets(AuxRing));  -- matrix in R
215
216QB;
217matrices := [TutEvalPolyInMatList(t, [A1,A2]) | t in QB];
218indent(matrices);
219C := BlockMat([[mat(AuxRing,M) * W | M in matrices]]);
220C;
221
222-- (4)
223D := det(C);  D;
224eval(D, [1,0,0,0]); --> not 0  --> OK
225--> [1,0,0,0] is a generator of QQ^3 as F-module
226V := ColMat([1,0,0,0]);
227
228C1000 := BlockMat([[M*V | M in matrices]]);
229invC := inverse(C1000);  invC;
230invC := mat(RelRing, invC);  invC;
231
232B1 := TutEvalPolyInMatList(ScalarProduct(QB, GetCol(invC,1)),
233			[A1,A2]);
234B2 := TutEvalPolyInMatList(ScalarProduct(QB, GetCol(invC,2)),
235			[A1,A2]);
236B3 := TutEvalPolyInMatList(ScalarProduct(QB, GetCol(invC,3)),
237			[A1,A2]);
238B4 := TutEvalPolyInMatList(ScalarProduct(QB, GetCol(invC,4)),
239			[A1,A2]);
240
241B1; B2; B3; B4;
242
243B1 * V;
244B2 * V;
245B3 * V;
246B4 * V;
247
248------------------------------------------------------------
249-- EXERCISE try with these matrices
250A1 := mat([[1,0,0],
251	   [0,1,0],
252	   [-1,-1,-1]]);
253A2 := mat([[0,1,0],
254	   [1,0,0],
255	   [0,0,1]]);
256------------------------------------------------------------
257-- UNIGENERATED FAMILIES
258
259A1 := matrix([[0,0,0,0,0,0,0],
260	      [0,0,0,0,0,-12,0],
261	      [0,0,0,0,0,0,-12],
262	      [1,0,0,0,0,0,0],
263	      [0,1,0,0,0,7,0],
264	      [0,0,0,1,0,0,0],
265	      [0,0,0,0,1,3,0]]);
266
267A2 := matrix([[0,0,0,0,0,0,0],
268	      [1,0,0,0,0,0,0],
269	      [0,1,0,0,0,0,0],
270	      [0,0,0,0,0,0,0],
271	      [0,0,0,1,0,0,0],
272	      [0,0,0,0,0,0,0],
273	      [0,0,0,0,0,1,0]]);
274A1*A2 = A2*A1;
275-- F commuting family generated by A1, A2
276
277use QQz;
278MinPoly(A1, z);  CharPoly(A1, z);  --> commendable
279MinPoly(A2, z);  CharPoly(A2, z);  --> non commendable
280
281-- Then now we know that F is unigenerated (by A1)
282--   in particular A2 is polynomial in A1:
283--   how do we find this polynomial?
284
285use P ::= QQ[x,y], lex;  --> lex is the key!
286Rel := BMForMat(P, [A2,A1]);  -- the generator last
287indent(Rel);
288
289f := x - Rel.GB[2];  f;
290TutEvalPolyInMatList(f, [A2,A1])  =  A2;
291------------------------------------------------------------
292