【技術分享】特徵值分解
- 2020 年 3 月 13 日
- 筆記
本文原作者：尹迪，經授權後發佈。
假設向量v是方陣A的特徵向量，可以表示成下面的形式：
這裡lambda表示特徵向量v所對應的特徵值。並且一個矩陣的一組特徵向量是一組正交向量。特徵值分解是將一個矩陣分解為下面的形式：
其中Q是這個矩陣A的特徵向量組成的矩陣。sigma是一個對角矩陣，每個對角線上的元素就是一個特徵值。
特徵值分解是一個提取矩陣特徵很不錯的方法，但是它只適合於方陣，對於非方陣，它不適合。這就需要用到奇異值分解。
1 源碼分析
MLlib使用ARPACK來求解特徵值分解。它的實現代碼如下
def symmetricEigs( mul: BDV[Double] => BDV[Double], n: Int, k: Int, tol: Double, maxIterations: Int): (BDV[Double], BDM[Double]) = { val arpack = ARPACK.getInstance() // tolerance used in stopping criterion val tolW = new doubleW(tol) // number of desired eigenvalues, 0 < nev < n val nev = new intW(k) // nev Lanczos vectors are generated in the first iteration // ncv-nev Lanczos vectors are generated in each subsequent iteration // ncv must be smaller than n val ncv = math.min(2 * k, n) // "I" for standard eigenvalue problem, "G" for generalized eigenvalue problem val bmat = "I" // "LM" : compute the NEV largest (in magnitude) eigenvalues val which = "LM" var iparam = new Array[Int](11) // use exact shift in each iteration iparam(0) = 1 // maximum number of Arnoldi update iterations, or the actual number of iterations on output iparam(2) = maxIterations // Mode 1: A*x = lambda*x, A symmetric iparam(6) = 1 var ido = new intW(0) var info = new intW(0) var resid = new Array[Double](n) var v = new Array[Double](n * ncv) var workd = new Array[Double](n * 3) var workl = new Array[Double](ncv * (ncv + 8)) var ipntr = new Array[Int](11) // call ARPACK's reverse communication, first iteration with ido = 0 arpack.dsaupd(ido, bmat, n, which, nev.`val`, tolW, resid, ncv, v, n, iparam, ipntr, workd, workl, workl.length, info) val w = BDV(workd) // ido = 99 : done flag in reverse communication while (ido.`val` != 99) { if (ido.`val` != -1 && ido.`val` != 1) { throw new IllegalStateException("ARPACK returns ido = " + ido.`val` + " This flag is not compatible with Mode 1: A*x = lambda*x, A symmetric.") } // multiply working vector with the matrix val inputOffset = ipntr(0) - 1 val outputOffset = ipntr(1) - 1 val x = w.slice(inputOffset, inputOffset + n) val y = w.slice(outputOffset, outputOffset + n) y := mul(x) // call ARPACK's reverse communication arpack.dsaupd(ido, bmat, n, which, nev.`val`, tolW, resid, ncv, v, n, iparam, ipntr, workd, workl, workl.length, info) } val d = new Array[Double](nev.`val`) val select = new Array[Boolean](ncv) // copy the Ritz vectors val z = java.util.Arrays.copyOfRange(v, 0, nev.`val` * n) // call ARPACK's post-processing for eigenvectors arpack.dseupd(true, "A", select, d, z, n, 0.0, bmat, n, which, nev, tol, resid, ncv, v, n, iparam, ipntr, workd, workl, workl.length, info) // number of computed eigenvalues, might be smaller than k val computed = iparam(4) val eigenPairs = java.util.Arrays.copyOfRange(d, 0, computed).zipWithIndex.map { r => (r._1, java.util.Arrays.copyOfRange(z, r._2 * n, r._2 * n + n)) } // sort the eigen-pairs in descending order val sortedEigenPairs = eigenPairs.sortBy(- _._1) // copy eigenvectors in descending order of eigenvalues val sortedU = BDM.zeros[Double](n, computed) sortedEigenPairs.zipWithIndex.foreach { r => val b = r._2 * n var i = 0 while (i < n) { sortedU.data(b + i) = r._1._2(i) i += 1 } } (BDV[Double](sortedEigenPairs.map(_._1)), sortedU) }
我們可以查看ARPACK的注釋詳細了解dsaupd和dseupd方法的作用。
