GL(n,p)矩阵的对角化 & Jordan Canonical Form

上上次说了一个模n矩阵的快速计算的两种方法，其实如果能把矩阵“对角化”的话，还可以进一步减少运算复杂度（不过这种方法在模n/模p下比较看运气）；另外，上次在说GL(n,p)的群阶时说了两个简略的证明，如果用Jordan Canonical Form来解释的话好像会更清晰.

普通矩阵的对角化·

这里说的“普通矩阵”，其实是指线性代数里学的$\mathbb{R}$上的矩阵，它的对角化方法在[GSLA]的6.2节有说到；给定一个$n*n$的非奇异矩阵$A$，要求它的对角矩阵$\Lambda$和对角化变化矩阵$X$的话：

1. 先求出$A$的特征多项式$p_A(x)$，根据Cayley–Hamilton定理：

   $$p_A(\lambda) = det(A-\lambda I_{n*n})$$

2. 通过解$p_A(\lambda)=0$可以解出$A$的$n$个特征值$\lambda_1,...,\lambda_n$，如果这$n$个特征值各不相同的话，就直接有：

   $$\Lambda = \begin{pmatrix} \lambda_1 & 0 & \cdots & 0 \\ 0 & \lambda_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda_n \\ \end{pmatrix}$$

3. 接着解出$n$个特征值对应的$n$个特征向量$x_1,...,x_n$，通过解方程：

   $$(A-\lambda_iI)x_i = 0$$

   可以解出各$x_i$，然后就有：

   $$X = [x_1, x_2, ..., x_n]$$

   可以简单验算一下有：

   $$AX = X\Lambda = [\lambda_1x_1, \lambda_2x_2, ..., \lambda_nx_n]$$

   最后对角化其实就是：

   $$A = X \Lambda X^{-1}$$

矩阵对角化的经典应用就是用来加快求幂：

$$A^k = X \Lambda^k X^{-1}$$

简要证明一下：原本$A^k$的话需要做$k$次（优化后$log(k)$次）$n*n$的矩阵乘法，每次复杂度是$O(n^3)$，（优化后是$O(n^{2+\alpha})$）；对角化后的$\Lambda^k$需要做$k$次（优化后$log(k)$次）$n*n$对角矩阵的乘法，每次复杂度$O(n)$，另外加两次矩阵乘法；所以在$k$很大的情况下是可以大大提高运算速度的：

$$O(log(k))*O(n)+2*O(n^{2+\alpha}) \approx O(n*log(k))$$

但是，在普通矩阵对角化的时候还有一些问题，比如：

1. 特征多项式$p_A(x)=0$无实数解的时候怎么办？
2. 解出的特征值有重解（$\lambda_i=\lambda_j$）的话怎么办？

域扩张（field extensions）·

对于第一个问题，可以用域扩张的方法，把实数域$\mathbb{R}$扩张到复数域$\mathbb{C}$，即原来$A$是在$\mathbb{R}$上做运算的，现在改成在$\mathbb{C}$上运算，那么$p_A(x)=0$就总会有解，只不过有的可能是复数解（有$i$），$\Lambda$和$X$上也会有复数。

若尔当型（Jordan Form）·

对于第二个问题，可以用Jordan Form来解决；Jordan Form（Jordan Matrix）就是$k$ 个块（Jordan Block）组成的块对角矩阵，如：

$$\begin{pmatrix} J_1 & 0 & \cdots & 0 \\ 0 & J_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & J_k \\ \end{pmatrix}$$

其中，$k$等于$A$的所有特征值的几何重数之和，也有说是独立的特征向量（independent eigenvectors）个数；一般情况下是$p_A(\lambda)$的不同根的个数，但也有相同的根$\lambda_i$对应多个Jordan Block的情况（，好像还需要看$(A-\lambda_iI)$对应零空间的维度？），可以参考下这里的一个例子；

Jordan Matrix中，$\lambda_i$对应的一个$J_i$长这样子：

$$\begin{pmatrix} \lambda_i & 1 & & & 0 \\ & \lambda_i & 1 & & \\ & & \ddots & \ddots \\ & & & \lambda_i & 1 \\ 0 & & & & \lambda_i \\ \end{pmatrix}$$

即对角线上是$\lambda_i$，对角线上一个（叫上对角线？）是$1$；其中每一个$J_i$的大小可以通过下面方式计算：

任何的非奇异方阵$A$都会相似于一个Jordan矩阵$J$，即可以找到一个$B$，使得：

$$A = BJB^{-1}$$

（顺便贴一下[GSLA] 8.3的说法，感觉比我说的简洁）

所以如果可以搞出$J$的话，也可以用类似对角化的方法加快幂运算：

$$A^k = B J^k B^{-1}$$

其中$J$做了$k$次方后仍然是只有“两条对角线”的的矩阵，即也是$k$次（优化后$log(k)$次）$O(n)$的矩阵乘法；

$J$的求法可以参考知乎的一个专栏；$B$的求法可以参考维基百科；大概就是：

1. 求$p_A(\lambda)$和各$\lambda_i$；

2. 求$\lambda_i$对应的Jordan块的阶数，根据阶数构造Jordan块，然后拼接Jordan矩阵；

3. 如果$J_i$的阶数为$1$的话，$B$中对应那一列就是$\lambda_i$对应的特征向量，即计算

   $$(A-\lambda_iI)b_i = 0$$

4. 如果$J_i$阶数为$l>1$的话，即可以算出$B$中的$l$个向量，其中第一个可以是$\lambda_i$的特征向量：

   $$(A-\lambda_iI)b_{i, 1} = 0$$

   然后其他剩余的$l-1$个向量可以通过解：

   $$(A-\lambda_iI)b_{i, j} = b_{i, j-1}$$

   得出（不过写代码时好像见到有无解的情况。。。）

GL(n, p)上的对角化·

$GL(n, \mathbb{F}_p)$上的对角化其实跟$\mathbb{R}$上的差不多，只不过有一些需要注意的问题。

求特征值·

按上面的说法，给定一个矩阵$A$，首先需要求它的特征多项式$p_A(\lambda)$，求的方法是一样的，即解：

$$p_A(\lambda) = det(A-\lambda I),\ \ in\ \mathbb{F}_p$$

只不过运算是在域$\mathbb{F}_p$上进行的；解特征值的方法也是一样的，即解：

$$p_A(\lambda)=0, \ \ in\ \mathbb{F}_p$$

然后这里就有了第一个问题：$\mathbb{F}_p$上的方程怎么解？

1. 首先可以对$p_A(\lambda)$进行分解，比如用[Ber70]的方法，但后来发现其实用Sage就可以直接分解，只不过要把域定义好：

   # Sage
   def factorFp(f, p):
     F.<a> = GF(p)
     R.<x> = PolynomialRing(F)
     return R(f).factor()

2. 如果$p_A(\lambda)$可以被分解为若干个度（degree）为$1$的多项式（即最高次数项是$\lambda$）的乘积，那么方程就可以直接被解了；比如在$\mathbb{F}_3$上有：

   $$p_A(\lambda)=\lambda^3+\lambda^2-\lambda-1=0, \ \ in\ \mathbb{F}_3$$

   可以被分解为：

   $$(\lambda+2) * (\lambda+1)^2 = 0, \ \ in\ \mathbb{F}_3$$

   就可以解出：

   $$\lambda = 1, 2, 2\ ; \ \ in\ \mathbb{F}_3$$

3. 如果$p_A(\lambda)$不能完全被分解为degree为$1$的多项式（或者直接不能被分解）的话，则需要使用域扩张的方法，扩展到更大的域上；比如在$\mathbb{F}_3$上有：

   $$p_A(\lambda)=\lambda^3+\lambda^2-\lambda+1=0, \ \ in\ \mathbb{F}_3$$

   在$\mathbb{F}_3$上$p_A(\lambda)$是不可约的，于是可以把$\mathbb{F}_3$扩张到：

   $$\mathbb{F}_{3^3}=\mathbb{F}_3[\lambda]/(\lambda^3+\lambda^2-\lambda+1)$$

   $\mathbb{F}_{3^3}$即模$p_A(\lambda)$的多项式环，其中系数的运算是在$\mathbb{F}_3$上的，记$\mathbb{F}_{3^3}$上的变量是$x$（只是个符号而已）；扩张后可以得到$p_A(\lambda)$的三个不同解（原理未知，咕）：

   $$\lambda=x^{3^0}, x^{3^1}, x^{3^2} ;\ \ in\ \mathbb{F}_{3^3}$$

   化简一下就是：

   $$\lambda = x, 2x^2+x+2, x^2+x ;\ \ in\ \mathbb{F}_{3^3}$$

4. 如果$p_A(\lambda)$被分解成多个degree大于1的不可约多项式的话，$\lambda$好像就会分布在多个不同的域上；比如在$\mathbb{F}_3$上有：

   $$p_A(\lambda)=\lambda^4+1=0, \ \ in\ \mathbb{F}_3$$

   $p_A(\lambda)$可以被分解为：

   $$p_A(\lambda) = (\lambda^2+\lambda+2) * (\lambda^2+2*\lambda+2), \ \ in\ \mathbb{F}_3$$

   $\lambda$就会分布在两个域上：$\mathbb{F}_3[\lambda]/(\lambda^2+\lambda+2)$和$\mathbb{F}_3[\lambda]/(\lambda^2+2*\lambda+2)$；

   感觉上应该会有方法扩张到同一个更大的域上（咕

对角化·

对角化的方法基本上是一样的，特征值$\lambda_i$各不相同的话就可以生成对角矩阵$\Lambda$，有重根的话就可以生成Jordan矩阵$J$；但像上面第4点说的，如果$\lambda$分布在不同的域上的话会有问题（咕

Demo Code：

# Sage
# 函数可以复制文章最下面的Demo Code

if __name__ == '__main__':
  p = 33285073849485750791903437807279991921
  m = matrix(IntegerModRing(p), [
    [1, 1, 1], 
    [1, 0, 0], 
    [0, 1, 0]])
  print(factorFp(m.charpoly(), p))
  B, J, Bi = jordanDiagon(m, p)
  print(p)
  print(J)
  print('--- ---')

  q = 38845988283830087557982578789883120419
  m = matrix(IntegerModRing(q), [
    [1, 1, 1], 
    [1, 0, 0], 
    [0, 1, 0]])
  print(factorFp(m.charpoly(), q))
  B, J, Bi = jordanDiagon(m, q)
  print(p)
  print(J)

之所以是Demo Code是因为上面问题是没有解决的，即不是每个矩阵都能分解，但维度$n \le 3$的话是都适用的；用了一个series中的矩阵和$e$分解的两个素数做例子，其中若$p=33285073849485750791903437807279991921$的话，可以解出：

$$\lambda_1 = 13232791622035946436448165355007395754, \ \ in\ \mathbb{F}_p \\ \lambda_2 = x, \ \ in\ \mathbb{F}_{p^2} \\ \lambda_3 = 33285073849485750791903437807279991920*y + 20052282227449804355455272452272596168, \ \ in\ \mathbb{F}_{p^2} \\$$

然后就可以构造出对角矩阵：

$$\Lambda= \begin{pmatrix} \lambda_1 & 0 & 0 \\ 0 & \lambda_2 & 0 \\ 0 & 0 & \lambda_3 \\ \end{pmatrix}$$

看上去同样是可以用这种方法，再求出$B$，然后用$B \Lambda^k B^{-1}$加快求幂，但是要注意扩展域后是在$\mathbb{F}_{p^2}$上做的乘法运算；在最坏的情况下，是在$\mathbb{F}_{p^n}$上（即$p_A(\lambda)$是不可约的情况）做$(n*log(k))$次乘法运算，一般的多项式乘法的复杂度是$O(n^2)$，使用FFT可优化到$O(n*log(n))$；即求$A$的$k$次幂可从原来的$O(n^{2+\alpha}*log(k))$优化到：

$$O(n^2*log(n)*log(k))$$

所以其实这种优化是比较“看脸”的，而且据说已经有矩阵乘法的方法可以做到$O(n^2*log(n))$了[Wied14]，不过如果$p_A(\lambda)$是可约的话还是可以优化一下的（

扩展·

既然现在知道$GL(n, \mathbb{F}_p)$的对角化方法了，那么就进一步解释上次计算$GL(n, \mathbb{F}_p)$的阶的证明了（大概讲一下成因，详细的证明不会写，咕）；另外这里也提供一种生成$GL(n, \mathbb{F}_p)$的“最大阶”（$p^n-1$）的方法：

证明的进一步解释·

首先说一下“最大阶”的证明，当时的说法是：对于$A \in GL(n, \mathbb{F}_p)$的话，最大的$A$的阶是$p^n-1$，原因是$p_A(\lambda)$最多$n$项，减去全$0$项的话最多就有$p^n-1$种特征多项式；当时论坛的那段回复后面还有一段：

大概意思就是$p_A(\lambda)$的解（在阶最大的情况？）是在扩张后的域$\mathbb{F}_{p^n}$上的，那么对角矩阵/Jordan矩阵的运算是在$\mathbb{F}_{p^n}$上的，先看一下比较简单的对角矩阵的情况：由于$\mathbb{F}_{p^n}$的阶就是$p^n-1$，所以$\Lambda$在做$p^n-1$次自乘后对角线上所有元素就会变为$1$，即：

$$\Lambda^{p^n-1} = I\ , \ \ \ \ in\ \mathbb{F}_{p^n}$$

Jordan矩阵的情况会比较复杂，经实验（懒得证了-），$J$的上对角线上有$1$的话，阶上就会产生因子$p$，在奇素数的情况下应该可以证到阶会比$p^n-1$小的？（注意如果$p_A(\lambda)$不可约的话特征值是$n$个$\mathbb{F}_{p^n}$上的不同值）， 咕

再说一下阶的证明，当时的说法是：

$$A^{p^r\prod_{i=1}^{n}\phi_i(p)} \equiv I \ \ (in\ GL_n(\mathbb{Z}_p))$$

先看$\phi_i(p)$的那$n$个因子，假设$p_A(\lambda)$可以被分解为$(n-i)$个degree为$1$的多项式和剩余的$1$个不可约多项式，那么$A$的特征值就有$(n-i)$个在$\mathbb{F}_{p}$上和$i$个在$\mathbb{F}_{p^i}$上，于是如果是特征值互不相等的情况的话，就有：

$$\Lambda^{p^i-1} = I\ , \ \ \ \ in\ \mathbb{F}_{p^i}$$

（注意$p^i-1$含因子$p-1$，也可看成$\mathbb{F}_{p}$可以被扩张到$\mathbb{F}_{p^i}$）；而$\phi_i(p)$其实是$p^i$的因子，所以就会有$\prod_{i=1}^{n}\phi_i(p)$这一块；

如果$A$的特征值互不相等的话，就是大概对角化到一个Jordan矩阵，会产生若干的$p$的因子，至于为什么是$p^r$可以看下[MW97]的第三章。。。

最大阶的生成元·

$A \in GL(n, \mathbb{F}_p)$的阶其实是和它的对角矩阵/Jordan矩阵的阶是一样的，按上面的说法，如果想达到最大阶的话，可以构造一个对角矩阵$\Lambda$，其元素都是$\mathbb{F}_{p^n}$上的，而且至少一个是$\mathbb{F}_{p^n}$的原根，那么$\Lambda$的阶就是$p^n-1$，即$A$的阶也是$p^n-1$；

根据上面的说法，$p_A(\lambda)$应该是一个$\mathbb{F}_p$上不可约的多项式（而且$p_A(\lambda)=0$至少有一个解是$\mathbb{F}_{p^n}$的原根），在假设能找到这样的$p_A(\lambda)$的前提下，可以通过以下方法找到$GL(n, \mathbb{F}_p)$中一个阶为$p^n-1$的元素：

1. 构造$p_A(\lambda)$的伴随矩阵，记为$C$：

   

2. 选取随机$R \in GL(n, \mathbb{F}_p)$，即可构造：

   $$A = RCR^{-1}$$

3. 根据伴随矩阵的性质，$C$的$n$个特征值就是$p_A(\lambda)$对应的$n$个解；根据上面说法，$C$的阶是$p^n-1$；$A≌C$即$A$的阶和$C$相同，也为$p^n-1$

Demo Code·

贴一个不太完善的参考代码：

# Sage
from collections import Counter

var('x')
def factorFp(f, p):
  F.<a> = GF(p)
  R.<x> = PolynomialRing(F)
  return R(f).factor()

def isIrreducible(f, p):
  ft = factorFp(f, p)
  return len(ft)==1 and ft[0][1]==1

# Actually, m.charpoly()
def getCharPoly(m):
  mm = m.change_ring(ZZ).change_ring(SR)  # Fp -> SR
  I = identity_matrix(mm.ncols())
  return (mm-x*I).det().expand()

def getLambdas(m, p):
  #cpm = getCharPoly(m)
  cpm = m.charpoly()
  fs = factorFp(cpm, p)
  F.<a> = GF(p)
  R.<x> = PolynomialRing(F)

  #lambdas = {}
  lambdas = {F: []}  # always F first ?
  try:
    for f in fs:
      for i in range(f[1]):
        deg = f[0].degree()
        if deg == 1:
          lam = F(-(f[0](x=0)))
          try:
            lambdas[F].append(lam)
          except:
            lambdas[F] = [lam]
        else:
          S = R.quotient(f[0], 'y')
          lambdas[S] = []
          lam = S(x)
          for d in range(deg):
            if d != 0:
              lam = power(lam, p)
            #lam = S(x^(p^d))
            lambdas[S].append(lam)
    return lambdas
  except Exception as e:
    # the right cpm will never get here ?
    print('???')
    print(e)
    exit(0)

# https://zhuanlan.zhihu.com/p/337122096
def jordanDiagon(m, p):  # default:right
  n = m.ncols()
  lambdas = getLambdas(m, p)

  if len(lambdas) > 2:
    raise NotImplementError("I don't know how to do that now ...")
  
  S = list(lambdas.keys())[-1]
  lams = []
  for Si in lambdas:
    lams += [S(x) for x in lambdas[Si]]
  lambdas = lams
  lambdas.sort()
  lams = Counter(lambdas)

  I = identity_matrix(S, n)
  Js = []
  evs = []  # eigenvectors
  for lam in lams.keys():
    numLam = lams[lam]
    m1 = m-lam*I
    numBlock = n - m1.rank()
    
    # get ri (rank is n for I)
    rs = [n]
    i = 1
    while True:
      ri = (m1^i).rank()
      if ri == rs[-1]:
        rs.append(ri)
        break
      rs.append(ri)
      i += 1

    # get orders
    jb = {}
    for i in range(len(rs))[1:-1]:
      tmp = rs[i-1] + rs[i+1] - 2*rs[i]
      if tmp != 0:
        jb[i] = tmp

    # there are k blocks with order d
    for di in range(len(jb.keys())):
      d = list(jb.keys())[di]
      for k in range(jb[d]):
        # for B
        ns = m1.right_kernel(basis='pivot')  # not work for reducible S ...
        if d==1:
          nb = ns.basis()
          evs.append(nb[di])
        else:  # TODO: no solution ???
          v0 = vector([0 for _ in range(n)])
          while True:
            ev0 = ns.random_element()
            if ev0==v0 or ev0 in evs:
              continue
            tmp = [ev0] 
            try:
              for i in range(1, d):
                tmp.append(m1 \ vector(tmp[-1]))
                #print(tmp)
            except:
              continue
            break
          evs += tmp
          
        # for J
        Ji = diagonal_matrix(S, [lam for _ in range(d)])
        jpos = 0
        for j in range(d):
          if j != d-1:
            Ji[jpos+j, jpos+j+1] = 1
        Js.append(Ji)
        jpos += d

  # Jordan block -> Jordan matrix
  #J = block_diagonal_matrix(S, Js)  # TypeError: Not a scalar type.
  J = zero_matrix(S, n)
  pos = 0
  for Ji in Js:
    ni = Ji.ncols()
    for ri in range(ni):
      for ci in range(ni):
        tmp = Ji[ri, ci]
        if tmp != 0:
          J[pos+ri, pos+ci] = tmp
    pos += ni

  B = matrix(S, evs).transpose()
  Binv = B.inverse()

  return B, J, Binv

def randMp(p, n):
  while True:
    m = matrix(IntegerModRing(p), [[randrange(p) for _ in range(n)] for _ in range(n)])
    if m.det()!=0:
      return m

# MA p194
# companion matrix of a polynomial (eigen values == roots of cp)
def companMp(cp, p):
  F.<a> = GF(p)
  R.<x> = PolynomialRing(F)
  f = R(cp)
  n = f.degree()
  m = zero_matrix(F, n)
  for i in range(1, n):
    m[i, i-1] = 1
  for i in range(n):
    m[i, n-1] = -f[i]
  return m

参考·

1. [GSLA] Introduction to Linear Algebra, Fifth Edition (2016), http://math.mit.edu/~gs/linearalgebra/
2. [Ber70] Berlekamp E R. Factoring polynomials over large finite fields[J]. Mathematics of computation, 1970, 24(111): 713-735.
3. [MW97] Menezes A J, Wu Y H. The discrete logarithm problem in GL (n, q)[J]. Ars Combinatoria, 1997, 47: 23-32.
4. [Wied14] Wiedermann J. Fast nondeterministic matrix multiplication via derandomization of Freivalds’ algorithm[C]//IFIP International Conference on Theoretical Computer Science. Springer, Berlin, Heidelberg, 2014: 123-135.
5. [MA] Horn R A, Johnson C R. Matrix analysis[M]. Cambridge university press, 2012.