【scala原理系列】 breeze.linalg.Breeze 特征值分解原理方法示例源码详解

原创已于 2023-12-26 13:11:07 修改 · 552 阅读

9 ·

CC 4.0 BY-SA版权

文章标签：

#scala

于 2023-12-22 15:15:35 首次发布

scala 专栏收录该内容

54 篇文章

订阅专栏

本文详细解释了Breeze.linalg包中Eig对象的原理，特别是Eig_DM_Impl隐式对象如何实现对DenseMatrix进行特征值分解，包括LAPACK库的调用和异常处理。并给出了示例代码展示了如何使用eig对象计算特征值和特征向量。

breeze.linalg.Breeze 特征值分解原理方法示例源码详解

文章目录

breeze.linalg.Breeze 特征值分解原理方法示例源码详解

原理

eig对象的源码实现了对任意矩阵进行特征值分解的功能。下面是关于eig对象的源码原理的详细解释：

1.`Eig` 类型定义：

case class Eig[V, M](eigenvalues: V, eigenvaluesComplex: V, eigenvectors: M)
type DenseEig = Eig[DenseVector[Double], DenseMatrix[Double]]

Eig 是一个包含特征值、复数特征值和特征向量的类，其中 V 表示特征值的类型，M 表示特征向量的类型。DenseEig 是一个具体的 Eig 类型，其中特征值和特征向量都是 DenseVector 和 DenseMatrix。

2.`Eig_DM_Impl` 隐式对象实现：

implicit object Eig_DM_Impl extends Impl[DenseMatrix[Double], DenseEig] {
  def apply(m: DenseMatrix[Double]): DenseEig = {
    // 检查矩阵是否符合要求

    // 分配空间

    // 调用 LAPACK 库进行特征值分解计算

    // 处理异常情况

    // 返回结果
  }
}

Eig_DM_Impl 是一个隐式对象，它实现了 Impl trait，并为 DenseMatrix[Double] 类型提供了特征值分解的功能。它的 apply 方法接受一个 DenseMatrix[Double] 类型的参数，并返回一个 DenseEig 类型的结果。

在 apply 方法中，首先对输入矩阵进行了一系列的检查，包括检查矩阵是否为空、是否为方阵以及是否存在 NaN 值等。然后，根据矩阵的尺寸分配了存储特征值和特征向量的空间。

接下来，调用 LAPACK 库中的 dgeev 函数进行特征值分解计算。这个函数会修改输入矩阵，并将结果存储在预先分配的空间中。

最后，处理了一些异常情况，如迭代次数超过限制或参数错误等。如果计算过程中出现问题，则抛出相应的异常。否则，将计算得到的特征值、复数特征值和特征向量作为 DenseEig 对象返回。

总之，eig 对象的源码实现了对任意矩阵进行特征值分解的功能。它使用 LAPACK 库中的函数进行计算，并通过一系列的检查和处理来确保计算的正确性和稳定性。

特征值分解（右特征向量）

此函数返回特征值的实部和虚部，以及相应的特征向量。对于大多数有趣的矩阵，
所有特征值的虚部都将为零（并且相应的特征向量将为实数）。任何复数特征值都会以共轭复数对的形式出现，
并且每个对应列的特征向量的实部和虚部将在特征向量矩阵中的相应列中。取复共轭以找到第二个特征向量。

示例

package org.example.scala

import breeze.linalg._
import breeze.numerics._

/**
 * 特征向量 特征值的求解
 */
object EigenvalueDecompositionDemo extends App {
  // 创建一个3x3的矩阵
  val matrix = DenseMatrix((1.0, 2.0, 3.0), (4.0, 5.0, 6.0), (7.0, 8.0, 9.0))

  // 使用eig对象计算特征值和特征向量
  val eigResult = eig(matrix)
  println("Eigenvalues:")
  println(eigResult.eigenvalues)
  println("Eigenvectors:")
  println(eigResult.eigenvectors)

  // 创建一个对称矩阵
  val symmetricMatrix = DenseMatrix((1.0, 2.0, 3.0), (2.0, 4.0, 5.0), (3.0, 5.0, 6.0))

  // 使用eigSym对象计算特征值和特征向量
  val eigSymResult = eigSym(symmetricMatrix)
  println("Eigenvalues:")
  println(eigSymResult.eigenvalues)
  println("Eigenvectors:")
  println(eigSymResult.eigenvectors)
}
//Eigenvalues:
//DenseVector(16.116843969807043, -1.1168439698070427, -1.3036777264747022E-15)
//Eigenvectors:
//-0.23197068724628617  -0.7858302387420671   0.40824829046386363
//-0.5253220933012336   -0.08675133925662833  -0.816496580927726
//-0.8186734993561815   0.61232756022881      0.4082482904638625
//Eigenvalues:
//DenseVector(-0.5157294715892579, 0.1709151888271789, 11.34481428276208)
//Eigenvectors:
//-0.7369762290995792  0.5910090485061029   -0.3279852776056817
//-0.3279852776056811  -0.7369762290995786  -0.5910090485061034
//0.5910090485061033   0.32798527760568236  -0.7369762290995783

用法

要使用eig对象进行特征值分解，可以按照以下步骤：

导入所需的包和类：

import breeze.linalg._

创建一个DenseMatrix[Double]类型的矩阵，表示要进行特征值分解的输入矩阵：

val matrix = DenseMatrix((1.0, 2.0, 3.0), (4.0, 5.0, 6.0), (7.0, 8.0, 9.0))

调用eig对象的apply方法，并传入输入矩阵作为参数，得到特征值和特征向量的结果：

val eigResult = eig(matrix)

可以通过eigResult对象访问计算得到的特征值和特征向量：

val eigenvalues = eigResult.eigenvalues
val eigenvectors = eigResult.eigenvectors

中文源码

/**
  * 特征值分解（右特征向量）
  *
  * 此函数返回特征值的实部和虚部，以及相应的特征向量。对于大多数有趣的矩阵，
  * 所有特征值的虚部都将为零（并且相应的特征向量将为实数）。任何复数特征值都会以共轭复数对的形式出现，
  * 并且每个对应列的特征向量的实部和虚部将在特征向量矩阵中的相应列中。取复共轭以找到第二个特征向量。
  *
  * 基于MTJ 0.9.12的EVD.java
  */
object eig extends UFunc {

  // TODO：可能我们应该只返回eigenValues：DV[Complex]？
  case class Eig[V, M](eigenvalues: V, eigenvaluesComplex: V, eigenvectors: M)
  type DenseEig = Eig[DenseVector[Double], DenseMatrix[Double]]


  implicit object Eig_DM_Impl extends Impl[DenseMatrix[Double], DenseEig] {
    def apply(m: DenseMatrix[Double]): DenseEig = {
      requireNonEmptyMatrix(m)
      requireSquareMatrix(m)
      require(!m.valuesIterator.exists(_.isNaN))

      val n = m.rows

      // 为分解分配空间
      val Wr = DenseVector.zeros[Double](n)
      val Wi = DenseVector.zeros[Double](n)

      val Vr = DenseMatrix.zeros[Double](n,n)

      // 寻找所需的工作空间
      val worksize = Array.ofDim[Double](1)
      val info = new intW(0)

      lapack.dgeev(
        "N", "V", n,
        Array.empty[Double], scala.math.max(1,n),
        Array.empty[Double], Array.empty[Double],
        Array.empty[Double], scala.math.max(1,n),
        Array.empty[Double], scala.math.max(1,n),
        worksize, -1, info)

      // 分配工作空间
      val lwork: Int = if (info.`val` != 0)
        scala.math.max(1,4*n)
      else
        scala.math.max(1,worksize(0).toInt)

      val work = Array.ofDim[Double](lwork)

      // 因子化！

      val A = DenseMatrix.zeros[Double](n, n)
      A := m
      lapack.dgeev(
        "N", "V", n,
        A.data, scala.math.max(1,n),
        Wr.data, Wi.data,
        Array.empty[Double], scala.math.max(1,n),
        Vr.data, scala.math.max(1,n),
        work, work.length, info)

      if (info.`val` > 0)
        throw new NotConvergedException(NotConvergedException.Iterations)
      else if (info.`val` < 0)
        throw new IllegalArgumentException()

      Eig(Wr, Wi, Vr)
    }
  }

}



/**
  * 计算给定实对称矩阵X的所有特征值（和可选的右特征向量）。
  */
object eigSym extends UFunc {
  case class EigSym[V, M](eigenvalues: V, eigenvectors: M)
  type DenseEigSym = EigSym[DenseVector[Double], DenseMatrix[Double]]
  implicit object EigSym_DM_Impl extends Impl[DenseMatrix[Double], DenseEigSym] {
    def apply(X: DenseMatrix[Double]): DenseEigSym = {
      doEigSym(X, true) match {
        case (ev, Some(rev)) => EigSym(ev, rev)
        case _ => throw new RuntimeException("Shouldn't be here!")
      }

    }
  }

  object justEigenvalues extends UFunc {
    implicit object EigSym_DM_Impl extends Impl[DenseMatrix[Double], DenseVector[Double]] {
      def apply(X: DenseMatrix[Double]): DenseVector[Double] = {
        doEigSym(X, false)._1
      }
    }

  }

  private def doEigSym(X: Matrix[Double], rightEigenvectors: Boolean): (DenseVector[Double], Option[DenseMatrix[Double]]) = {
    requireNonEmptyMatrix(X)

    // 由于LAPACK不检查给定矩阵是否对称，因此我们必须在此处进行检查（或者去掉这个时间浪费，
    // 只要调用此函数的调用方清楚地知道只使用给定矩阵的下三角部分，并且没有检查对称性）。
    requireSymmetricMatrix(X)

    // 复制X的下三角部分。LAPACK将结果存储在A中。
    val A     = lowerTriangular(X)

    val N     = X.rows
    val evs   = DenseVector.zeros[Double](N)
    val lwork = scala.math.max(1, 3*N-1)
    val work  = Array.ofDim[Double](lwork)
    val info  = new intW(0)
    lapack.dsyev(
      if (rightEigenvectors) "V" else "N" /* eigenvalues N, eigenvalues & eigenvectors "V" */,
      "L" /* lower triangular */,
      N /* number of rows */, A.data, scala.math.max(1, N) /* LDA */,
      evs.data,
      work /* workspace */, lwork /* workspace size */,
      info
    )
    // 如果info.`val` < 0，则告诉我们调用dsyev的第i个参数错误（其中i == |info.`val`|）。
    assert(info.`val` >= 0)

    if (info.`val` > 0)
      throw new NotConvergedException(NotConvergedException.Iterations)

    (evs, if (rightEigenvectors) Some(A) else None)
  }
}