特征值分解

  假设向量v是方阵A的特征向量，可以表示成下面的形式：

1.1

  这里lambda表示特征向量v所对应的特征值。并且一个矩阵的一组特征向量是一组正交向量。特征值分解是将一个矩阵分解为下面的形式：

1.2

  其中Q是这个矩阵A的特征向量组成的矩阵。sigma是一个对角矩阵，每个对角线上的元素就是一个特征值。

  特征值分解是一个提取矩阵特征很不错的方法，但是它只适合于方阵，对于非方阵，它不适合。这就需要用到奇异值分解。

1 源码分析

  MLlib使用ARPACK来求解特征值分解。它的实现代码如下

  1. def symmetricEigs(
  2.       mul: BDV[Double] => BDV[Double],
  3.       n: Int,
  4.       k: Int,
  5.       tol: Double,
  6.       maxIterations: Int): (BDV[Double], BDM[Double]) = {
  7.     val arpack = ARPACK.getInstance()
  8.     // tolerance used in stopping criterion
  9.     val tolW = new doubleW(tol)
  10.     // number of desired eigenvalues, 0 < nev < n
  11.     val nev = new intW(k)
  12.     // nev Lanczos vectors are generated in the first iteration
  13.     // ncv-nev Lanczos vectors are generated in each subsequent iteration
  14.     // ncv must be smaller than n
  15.     val ncv = math.min(2 * k, n)
  16.     // "I" for standard eigenvalue problem, "G" for generalized eigenvalue problem
  17.     val bmat = "I"
  18.     // "LM" : compute the NEV largest (in magnitude) eigenvalues
  19.     val which = "LM"
  20.     var iparam = new Array[Int](11)
  21.     // use exact shift in each iteration
  22.     iparam(0) = 1
  23.     // maximum number of Arnoldi update iterations, or the actual number of iterations on output
  24.     iparam(2) = maxIterations
  25.     // Mode 1: A*x = lambda*x, A symmetric
  26.     iparam(6) = 1
  28.     var ido = new intW(0)
  29.     var info = new intW(0)
  30.     var resid = new Array[Double](n)
  31.     var v = new Array[Double](n * ncv)
  32.     var workd = new Array[Double](n * 3)
  33.     var workl = new Array[Double](ncv * (ncv + 8))
  34.     var ipntr = new Array[Int](11)
  36.     // call ARPACK's reverse communication, first iteration with ido = 0
  37.     arpack.dsaupd(ido, bmat, n, which, nev.`val`, tolW, resid, ncv, v, n, iparam, ipntr, workd,
  38.       workl, workl.length, info)
  39.     val w = BDV(workd)
  40.     // ido = 99 : done flag in reverse communication
  41.     while (ido.`val` != 99) {
  42.       if (ido.`val` != -1 && ido.`val` != 1) {
  43.         throw new IllegalStateException("ARPACK returns ido = " + ido.`val` +
  44.             " This flag is not compatible with Mode 1: A*x = lambda*x, A symmetric.")
  45.       }
  46.       // multiply working vector with the matrix
  47.       val inputOffset = ipntr(0) - 1
  48.       val outputOffset = ipntr(1) - 1
  49.       val x = w.slice(inputOffset, inputOffset + n)
  50.       val y = w.slice(outputOffset, outputOffset + n)
  51.       y := mul(x)
  52.       // call ARPACK's reverse communication
  53.       arpack.dsaupd(ido, bmat, n, which, nev.`val`, tolW, resid, ncv, v, n, iparam, ipntr,
  54.         workd, workl, workl.length, info)
  55.     }
  57.     val d = new Array[Double](nev.`val`)
  58.     val select = new Array[Boolean](ncv)
  59.     // copy the Ritz vectors
  60.     val z = java.util.Arrays.copyOfRange(v, 0, nev.`val` * n)
  62.     // call ARPACK's post-processing for eigenvectors
  63.     arpack.dseupd(true, "A", select, d, z, n, 0.0, bmat, n, which, nev, tol, resid, ncv, v, n,
  64.       iparam, ipntr, workd, workl, workl.length, info)
  66.     // number of computed eigenvalues, might be smaller than k
  67.     val computed = iparam(4)
  69.     val eigenPairs = java.util.Arrays.copyOfRange(d, 0, computed).zipWithIndex.map { r =>
  70.       (r._1, java.util.Arrays.copyOfRange(z, r._2 * n, r._2 * n + n))
  71.     }
  73.     // sort the eigen-pairs in descending order
  74.     val sortedEigenPairs = eigenPairs.sortBy(- _._1)
  76.     // copy eigenvectors in descending order of eigenvalues
  77.     val sortedU = BDM.zeros[Double](n, computed)
  78.     sortedEigenPairs.zipWithIndex.foreach { r =>
  79.       val b = r._2 * n
  80.       var i = 0
  81.       while (i < n) {
  82.         sortedU.data(b + i) = r._1._2(i)
  83.         i += 1
  84.       }
  85.     }
  86.     (BDV[Double](sortedEigenPairs.map(_._1)), sortedU)
  87.   }

  我们可以查看ARPACK的注释详细了解dsaupd和dseupd方法的作用。