Improve Jacobi algorithm readability by extracting some logic into helper fun
This commit is contained in:
parent
7a72a0b979
commit
7aff774bc1
@ -10,6 +10,7 @@ package space.kscience.kmath.tensors.core
|
|||||||
|
|
||||||
import space.kscience.kmath.misc.PerformancePitfall
|
import space.kscience.kmath.misc.PerformancePitfall
|
||||||
import space.kscience.kmath.nd.MutableStructure2D
|
import space.kscience.kmath.nd.MutableStructure2D
|
||||||
|
import space.kscience.kmath.nd.Structure2D
|
||||||
import space.kscience.kmath.nd.StructureND
|
import space.kscience.kmath.nd.StructureND
|
||||||
import space.kscience.kmath.nd.as1D
|
import space.kscience.kmath.nd.as1D
|
||||||
import space.kscience.kmath.nd.as2D
|
import space.kscience.kmath.nd.as2D
|
||||||
@ -885,7 +886,7 @@ public open class DoubleTensorAlgebra :
|
|||||||
return Triple(uTensor.transpose(), sTensor, vTensor.transpose())
|
return Triple(uTensor.transpose(), sTensor, vTensor.transpose())
|
||||||
}
|
}
|
||||||
|
|
||||||
override fun StructureND<Double>.symEig(): Pair<DoubleTensor, DoubleTensor> = symEig(epsilon = 1e-15)
|
override fun StructureND<Double>.symEig(): Pair<DoubleTensor, DoubleTensor> = symEigJacobi(epsilon = 1e-10)
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Returns eigenvalues and eigenvectors of a real symmetric matrix input or a batch of real symmetric matrices,
|
* Returns eigenvalues and eigenvectors of a real symmetric matrix input or a batch of real symmetric matrices,
|
||||||
@ -922,96 +923,99 @@ public open class DoubleTensorAlgebra :
|
|||||||
return eig to v
|
return eig to v
|
||||||
}
|
}
|
||||||
|
|
||||||
// TODO
|
public fun StructureND<Double>.symEigJacobi(epsilon: Double): Pair<DoubleTensor, DoubleTensor> {
|
||||||
// 1. Cyclic method
|
|
||||||
// 2. Sort eigenvalues
|
|
||||||
public fun StructureND<Double>.symEig(epsilon: Double): Pair<DoubleTensor, DoubleTensor> {
|
|
||||||
checkSymmetric(tensor, epsilon)
|
checkSymmetric(tensor, epsilon)
|
||||||
|
|
||||||
val ii = tensor.minusIndex(-2)
|
|
||||||
val jj = tensor.minusIndex(-1)
|
|
||||||
val n = tensor.numElements
|
|
||||||
|
|
||||||
val size = this.dimension
|
val size = this.dimension
|
||||||
val commonShape = this.shape.sliceArray(0 until size - 2) + intArrayOf(1, 1)
|
val s = zeros(this.shape)
|
||||||
|
val eigenvalues = zeros(this.shape.sliceArray(0 until size - 1))
|
||||||
|
|
||||||
var d = this.copy()
|
var eigenvalueStart = 0
|
||||||
var s = diagonalEmbedding(ones(this.shape.sliceArray(0 until size - 1)))
|
var eigenvectorStart = 0
|
||||||
|
for (matrix in tensor.matrixSequence()) {
|
||||||
|
matrix.as2D().jacobiHelper(eigenvalues, s, eigenvalueStart, eigenvectorStart, epsilon)
|
||||||
|
eigenvalueStart += this.shape.last()
|
||||||
|
eigenvectorStart += this.shape.last() * this.shape.last()
|
||||||
|
}
|
||||||
|
|
||||||
|
// TODO sort eigenvalues
|
||||||
|
return eigenvalues to s
|
||||||
|
}
|
||||||
|
|
||||||
|
private fun MutableStructure2D<Double>.jacobiHelper(
|
||||||
|
eigenvalues: DoubleTensor,
|
||||||
|
eigenvectors: DoubleTensor,
|
||||||
|
eigenvalueStart: Int,
|
||||||
|
eigenvectorStart: Int,
|
||||||
|
epsilon: Double
|
||||||
|
) {
|
||||||
|
var d = this
|
||||||
|
var s = eye(this.shape[0])
|
||||||
|
|
||||||
|
// TODO implement cyclic method
|
||||||
do {
|
do {
|
||||||
// 1. Find max element by abs value that is not on diagonal
|
// 1. Find max element by abs value that is not on diagonal
|
||||||
val buffer = MutableBuffer.boxing(commonShape.reduce(Int::times)) { Triple(0.0, 0, 0) }
|
var maxOffDiagonalElement = 0.0
|
||||||
val maxOffDiagonalElements = BufferedTensor(commonShape, buffer, 0)
|
var maxElementIndex = Pair(0, 0)
|
||||||
for (offset in 0 until n) {
|
for (i in 0 until this.rowNum) {
|
||||||
val multiIndex = d.linearStructure.index(offset)
|
for (j in 0 until this.colNum) {
|
||||||
if (multiIndex[ii] != multiIndex[jj]) {
|
if (i == j) continue
|
||||||
val value = d.mutableBuffer.array()[offset]
|
if (abs(d[i, j]) > maxOffDiagonalElement) {
|
||||||
|
maxOffDiagonalElement = abs(d[i, j])
|
||||||
val commonIndex = multiIndex.sliceArray(0 until size - 2) + intArrayOf(0, 0)
|
maxElementIndex = i to j
|
||||||
if (abs(value) > maxOffDiagonalElements[commonIndex].first) {
|
|
||||||
maxOffDiagonalElements[commonIndex] = Triple(abs(value), multiIndex[ii], multiIndex[jj])
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
// 2. Evaluate "rotation" angle
|
// 2. Evaluate "rotation" angle
|
||||||
val angles = zeros(commonShape)
|
val dIJ = d[maxElementIndex.first, maxElementIndex.second]
|
||||||
for (offset in 0 until maxOffDiagonalElements.numElements) {
|
val dII = d[maxElementIndex.first, maxElementIndex.first]
|
||||||
val (_, i, j) = maxOffDiagonalElements.mutableBuffer[offset]
|
val dJJ = d[maxElementIndex.second, maxElementIndex.second]
|
||||||
val multiIndex = maxOffDiagonalElements.linearStructure.index(offset)
|
|
||||||
|
|
||||||
val dIJ = d.mutableBuffer[d.linearStructure.offset(multiIndex.also { it[ii] = i; it[jj] = j })]
|
val angle = if (dII == dJJ) {
|
||||||
val dII = d.mutableBuffer[d.linearStructure.offset(multiIndex.also { it[ii] = i; it[jj] = i })]
|
if (dIJ > 0) PI / 4 else -PI / 4
|
||||||
val dJJ = d.mutableBuffer[d.linearStructure.offset(multiIndex.also { it[ii] = j; it[jj] = j })]
|
} else {
|
||||||
|
0.5 * atan(2 * dIJ / (dJJ - dII))
|
||||||
angles.mutableBuffer.array()[offset] = if (dII == dJJ) {
|
|
||||||
if (dIJ > 0) PI / 4 else -PI / 4
|
|
||||||
} else {
|
|
||||||
0.5 * atan(2 * dIJ / (dJJ - dII))
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|
||||||
// 3. Build rotation tensor `s1`
|
// 3. Build rotation tensor `s1`
|
||||||
val s1 = diagonalEmbedding(ones(this.shape.sliceArray(0 until size - 1)))
|
val s1 = eye(this.rowNum)
|
||||||
for (offset in 0 until n) {
|
for (i in 0 until this.rowNum) {
|
||||||
val multiIndex = d.linearStructure.index(offset)
|
for (j in 0 until this.colNum) {
|
||||||
|
s1.mutableBuffer.array()[i * this.rowNum + j] = when {
|
||||||
val commonIndex = multiIndex.sliceArray(0 until size - 2) + intArrayOf(0, 0)
|
maxElementIndex.first == i && maxElementIndex.first == j -> cos(angle)
|
||||||
val (_, i, j) = maxOffDiagonalElements[commonIndex]
|
maxElementIndex.second == i && maxElementIndex.second == j -> cos(angle)
|
||||||
val angleValue = angles[commonIndex]
|
maxElementIndex.first == i && maxElementIndex.second == j -> sin(angle)
|
||||||
s1.mutableBuffer.array()[offset] = when {
|
maxElementIndex.first == j && maxElementIndex.second == i -> -sin(angle)
|
||||||
multiIndex[ii] == i && multiIndex[jj] == i || multiIndex[ii] == j && multiIndex[jj] == j -> cos(angleValue)
|
else -> s1.mutableBuffer.array()[i * this.rowNum + j]
|
||||||
multiIndex[ii] == i && multiIndex[jj] == j -> sin(angleValue)
|
}
|
||||||
multiIndex[ii] == j && multiIndex[jj] == i -> -sin(angleValue)
|
|
||||||
else -> s1.mutableBuffer.array()[offset]
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
// 4. Evaluate new tensor
|
// 4. Evaluate new tensor
|
||||||
d = (s1.transpose() dot d) dot s1
|
d = ((s1.transpose() dot d) dot s1).as2D()
|
||||||
s = s dot s1
|
s = s dot s1
|
||||||
if (d.isDiagonal(epsilon)) break
|
if (d.isDiagonal(epsilon)) break
|
||||||
} while(true)
|
} while(true)
|
||||||
|
|
||||||
val eigenvalues = zeros(d.shape.sliceArray(0 until size - 1))
|
// 5. Copy result
|
||||||
for (offset in 0 until n) {
|
for (i in 0 until this.rowNum) {
|
||||||
val multiIndex = d.linearStructure.index(offset)
|
for (j in 0 until this.colNum) {
|
||||||
if (multiIndex[ii] == multiIndex[jj]) {
|
eigenvectors.mutableBuffer.array()[eigenvectorStart + i * this.rowNum + j] = s.mutableBuffer.array()[i * this.rowNum + j]
|
||||||
eigenvalues[multiIndex.sliceArray(0 until size - 1)] = d.mutableBuffer.array()[offset]
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
return eigenvalues to s
|
for (i in 0 until this.rowNum) {
|
||||||
|
eigenvalues.mutableBuffer.array()[eigenvalueStart + i] = d[i, i]
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
public fun StructureND<Double>.isDiagonal(epsilon: Double = 1e-9): Boolean {
|
public fun Structure2D<Double>.isDiagonal(epsilon: Double = 1e-9): Boolean {
|
||||||
val ii = tensor.minusIndex(-2)
|
for (i in 0 until this.rowNum) {
|
||||||
val jj = tensor.minusIndex(-1)
|
for (j in 0 until this.colNum) {
|
||||||
|
if (i != j && abs(this[i, j]) > epsilon) {
|
||||||
for (offset in 0 until tensor.numElements) {
|
return false
|
||||||
val multiIndex = tensor.linearStructure.index(offset)
|
}
|
||||||
if (multiIndex[ii] != multiIndex[jj] && abs(tensor.mutableBuffer.array()[offset]) > epsilon) {
|
|
||||||
return false
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
Loading…
Reference in New Issue
Block a user