雅克比矩阵的scala实现

在向量分析中, 雅可比矩阵是一阶偏导数以一定方式排列成的矩阵, 其行列式称为雅可比行列式. 还有, 在代数几何中, 代数曲线的雅可比量表示雅可比簇：伴随该曲线的一个代数群, 曲线可以嵌入其中。

矩阵的特征值和特征向量是线性代数以及矩阵论中非常重要的一个概念。在遥感领域也是经常用到，比如多光谱以及高光谱图像的主成分分析要求解波段间协方差矩阵或者相关系数矩阵的特征值和特征向量。

根据普通线性代数中的概念，特征值和特征向量可以用传统的方法求得，但是实际项目中一般都是用数值分析的方法来计算。

雅克比方法用于求实对称阵的全部特征值、特征向量。对于实对称阵 A，必有正交阵 U，使U TA U = D。其中 D 是对角阵，其主对角线元 li 是 A 的特征值. 正交阵 U 的第 j 列是 A 的属于 li 的特征向量。Jacobi 方法用平面旋转对矩阵 A 做相似变换，化A 为对角阵，进而求出特征值与特征向量。

关于雅克比矩阵计算的原理不多做介绍，具体可参考博客：http://blog.csdn.net/zhouxuguang236/article/details/40212143

具体的scala实现如下：

def jacbicor(accuracy:Double,iternum:Int) = {
val matrixVector:Array[Array[Double]] = Array.ofDim[Double](rownum,colnum)
for(i <- 0 until rownum){
for(j <- 0 until colnum){
if(i == j){
matrixVector(i)(j) = 1.0
}else{
matrixVector(i)(j) = 0.0
}
}
}
var nCount = 0
var rotation = true
while(rotation){
var maxMatrixElement = this.matrix(0)(1)
var nRow = 0
var nCol = 1
for(i <- 0 until this.rownum){
for( j <- 0 until this.colnum){
var matrixElement = Math.abs(this.matrix(i)(j))
if((i != j)&&(matrixElement > maxMatrixElement)){
maxMatrixElement = matrixElement
nRow = i
nCol = j
}
}
}
if((nCount == iternum)||(maxMatrixElement < accuracy)){
rotation = false
}
nCount += 1
val matrixApp = this.matrix(nRow)(nRow)
val matrixApq = this.matrix(nRow)(nCol)
val matrixAqq = this.matrix(nCol)(nCol)
//count Rotation
val matrixAngle = 0.5*Math.atan2(-2*matrixApq,matrixAqq - matrixApp)
val matrixSinTheta = Math.sin(matrixAngle)
val matrixCosTheta = Math.cos(matrixAngle)
val matrixSin2Theta = Math.sin(2*matrixAngle)
val matrixCos2Theta = Math.cos(2*matrixAngle)
//matrix iter
this.matrix(nRow)(nRow) = matrixApp*matrixCosTheta*matrixCosTheta +
matrixAqq*matrixSinTheta*matrixSinTheta + 2*matrixApq*matrixCosTheta*matrixSinTheta
this.matrix(nCol)(nCol) = matrixApp*matrixSinTheta*matrixSinTheta +
matrixAqq*matrixCosTheta*matrixCosTheta - 2*matrixApq*matrixCosTheta*matrixSinTheta
this.matrix(nRow)(nCol) = 0.5*(matrixAqq-matrixApp)*matrixSin2Theta + matrixApq*matrixCos2Theta
this.matrix(nCol)(nRow) = this.matrix(nRow)(nCol)
for(i <- 0 until this.rownum){
if((i != nRow)&&( i != nCol)){
maxMatrixElement = this.matrix(i)(nRow)
this.matrix(i)(nRow) = this.matrix(i)(nCol)*matrixSinTheta + maxMatrixElement*matrixCosTheta
this.matrix(i)(nCol) = this.matrix(i)(nCol)*matrixCosTheta - maxMatrixElement*matrixSinTheta
}
}
for(j <- 0 until this.colnum){
if((j != nRow)&&(j != nCol)){
maxMatrixElement = this.matrix(nRow)(j)
this.matrix(nRow)(j) = this.matrix(nCol)(j)*matrixSinTheta + maxMatrixElement*matrixCosTheta
this.matrix(nCol)(j) = this.matrix(nCol)(j)*matrixCosTheta - maxMatrixElement*matrixSinTheta
}
}
for(i <- 0 until this.rownum){
maxMatrixElement = this.matrix(i)(nRow)
matrixVector(i)(nRow) = matrixVector(i)(nCol)*matrixSinTheta + maxMatrixElement*matrixCosTheta
matrixVector(i)(nCol) = matrixVector(i)(nCol)*matrixCosTheta - maxMatrixElement*matrixSinTheta
}
} //while end
for(j <- 0 until this.colnum){
var sumVector = 0.0
for(i <- 0 until this.rownum){
sumVector += matrixVector(i)(j)
}
if(sumVector < 0){
for(i <- 0 until this.colnum){
matrixVector(i)(j) *= -1
}
}
}
var eigvalueVectorMap = scala.collection.mutable.Map[Double,Array[Double]]()
for(j <- 0 until this.colnum){
val eigValue = this.matrix(j)(j)
val eigVector = ArrayBuffer[Double]()
for(i <- 0 until this.rownum){
eigVector += matrixVector(i)(j)
}
eigvalueVectorMap += (eigValue -> eigVector.toArray)
}
eigvalueVectorMap
}