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)
}
}