降维 - EVD(特征值分解)

优质
小牛编辑
128浏览
2023-12-01

  假设向量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
  27. var ido = new intW(0)
  28. var info = new intW(0)
  29. var resid = new Array[Double](n)
  30. var v = new Array[Double](n * ncv)
  31. var workd = new Array[Double](n * 3)
  32. var workl = new Array[Double](ncv * (ncv + 8))
  33. var ipntr = new Array[Int](11)
  34. // call ARPACK's reverse communication, first iteration with ido = 0
  35. arpack.dsaupd(ido, bmat, n, which, nev.`val`, tolW, resid, ncv, v, n, iparam, ipntr, workd,
  36. workl, workl.length, info)
  37. val w = BDV(workd)
  38. // ido = 99 : done flag in reverse communication
  39. while (ido.`val` != 99) {
  40. if (ido.`val` != -1 && ido.`val` != 1) {
  41. throw new IllegalStateException("ARPACK returns ido = " + ido.`val` +
  42. " This flag is not compatible with Mode 1: A*x = lambda*x, A symmetric.")
  43. }
  44. // multiply working vector with the matrix
  45. val inputOffset = ipntr(0) - 1
  46. val outputOffset = ipntr(1) - 1
  47. val x = w.slice(inputOffset, inputOffset + n)
  48. val y = w.slice(outputOffset, outputOffset + n)
  49. y := mul(x)
  50. // call ARPACK's reverse communication
  51. arpack.dsaupd(ido, bmat, n, which, nev.`val`, tolW, resid, ncv, v, n, iparam, ipntr,
  52. workd, workl, workl.length, info)
  53. }
  54. val d = new Array[Double](nev.`val`)
  55. val select = new Array[Boolean](ncv)
  56. // copy the Ritz vectors
  57. val z = java.util.Arrays.copyOfRange(v, 0, nev.`val` * n)
  58. // call ARPACK's post-processing for eigenvectors
  59. arpack.dseupd(true, "A", select, d, z, n, 0.0, bmat, n, which, nev, tol, resid, ncv, v, n,
  60. iparam, ipntr, workd, workl, workl.length, info)
  61. // number of computed eigenvalues, might be smaller than k
  62. val computed = iparam(4)
  63. val eigenPairs = java.util.Arrays.copyOfRange(d, 0, computed).zipWithIndex.map { r =>
  64. (r._1, java.util.Arrays.copyOfRange(z, r._2 * n, r._2 * n + n))
  65. }
  66. // sort the eigen-pairs in descending order
  67. val sortedEigenPairs = eigenPairs.sortBy(- _._1)
  68. // copy eigenvectors in descending order of eigenvalues
  69. val sortedU = BDM.zeros[Double](n, computed)
  70. sortedEigenPairs.zipWithIndex.foreach { r =>
  71. val b = r._2 * n
  72. var i = 0
  73. while (i < n) {
  74. sortedU.data(b + i) = r._1._2(i)
  75. i += 1
  76. }
  77. }
  78. (BDV[Double](sortedEigenPairs.map(_._1)), sortedU)
  79. }

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