kscience/kmath
9
0
Fork 1
Code Issues Pull Requests 6 Packages Projects Releases 16 Wiki Activity

levenbergMarquardt cleanup

Browse Source
This commit is contained in:
Alexander Nozik 2023-07-28 20:56:31 +03:00
parent 1e2a8a40e5
commit 976714475e
2 changed files with 31 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/LevenbergMarquardtAlgorithm.kt
View File

@ -120,8 +120,6 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI
0.0, inputData.startParameters, TypeOfConvergence.NoConvergence 0.0, inputData.startParameters, TypeOfConvergence.NoConvergence
) )
val eps = 2.2204e-16
val settings = LMSettings(0, 0, inputData.exampleNumber) val settings = LMSettings(0, 0, inputData.exampleNumber)
settings.funcCalls = 0 // running count of function evaluations settings.funcCalls = 0 // running count of function evaluations
@ -214,7 +212,7 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI
stop = true stop = true
} }
var lambda = 1.0 var lambda: Double
var nu = 1 var nu = 1
if (updateType == 1) { if (updateType == 1) {
@ -273,8 +271,8 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI
val alpha = 1.0 val alpha = 1.0
if (updateType == 2) { // Quadratic if (updateType == 2) { // Quadratic
// One step of quadratic line update in the h direction for minimum X2 // One step of quadratic line update in the h direction for minimum X2
val alpha = (jtWdy.transpose() dot h) / ((X2Try - x2) / 2.0 + 2 * (jtWdy.transpose() dot h)) val alphaTensor = (jtWdy.transpose() dot h) / ((X2Try - x2) / 2.0 + 2 * (jtWdy.transpose() dot h))
h = h dot alpha h = h dot alphaTensor
pTry = (p + h).as2D() // update only [idx] elements pTry = (p + h).as2D() // update only [idx] elements
pTry = smallestElementComparison( pTry = smallestElementComparison(
largestElementComparison(minParameters, pTry), largestElementComparison(minParameters, pTry),
@ -388,7 +386,7 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI
resultInfo.typeOfConvergence = TypeOfConvergence.InGradient resultInfo.typeOfConvergence = TypeOfConvergence.InGradient
stop = true stop = true
} }
if ((abs(h.as2D()).div(abs(p) + 1e-12)).max() < epsilon2 && settings.iteration > 2) { if ((abs(h.as2D()) / (abs(p) + 1e-12)).max() < epsilon2 && settings.iteration > 2) {
resultInfo.typeOfConvergence = TypeOfConvergence.InParameters resultInfo.typeOfConvergence = TypeOfConvergence.InParameters
stop = true stop = true
} }

49
kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt
View File

@ -7,7 +7,10 @@ package space.kscience.kmath.tensors.core.internal
import space.kscience.kmath.nd.* import space.kscience.kmath.nd.*
import space.kscience.kmath.operations.invoke import space.kscience.kmath.operations.invoke
import space.kscience.kmath.structures.* import space.kscience.kmath.structures.DoubleBuffer
import space.kscience.kmath.structures.IntBuffer
import space.kscience.kmath.structures.asBuffer
import space.kscience.kmath.structures.indices
import space.kscience.kmath.tensors.core.* import space.kscience.kmath.tensors.core.*
import kotlin.math.abs import kotlin.math.abs
import kotlin.math.max import kotlin.math.max
@ -329,14 +332,16 @@ private fun SIGN(a: Double, b: Double): Double {
return -abs(a) return -abs(a)
} }
internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2D<Double>, w: BufferedTensor<Double>, internal fun MutableStructure2D<Double>.svdGolubKahanHelper(
v: MutableStructure2D<Double>, iterations: Int, epsilon: Double) { u: MutableStructure2D<Double>, w: BufferedTensor<Double>,
v: MutableStructure2D<Double>, iterations: Int, epsilon: Double,
) {
val shape = this.shape val shape = this.shape
val m = shape.component1() val m = shape.component1()
val n = shape.component2() val n = shape.component2()
var f = 0.0 var f: Double
val rv1 = DoubleArray(n) val rv1 = DoubleArray(n)
var s = 0.0 var s: Double
var scale = 0.0 var scale = 0.0
var anorm = 0.0 var anorm = 0.0
var g = 0.0 var g = 0.0
@ -362,10 +367,10 @@ internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2
s += this[k, i] * this[k, i] s += this[k, i] * this[k, i]
} }
f = this[i, i] f = this[i, i]
if (f >= 0) { g = if (f >= 0) {
g = (-1) * abs(sqrt(s)) -abs(sqrt(s))
} else { } else {
g = abs(sqrt(s)) abs(sqrt(s))
} }
val h = f * g - s val h = f * g - s
this[i, i] = f - g this[i, i] = f - g
@ -402,10 +407,10 @@ internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2
s += this[i, k] * this[i, k] s += this[i, k] * this[i, k]
} }
f = this[i, l] f = this[i, l]
if (f >= 0) { g = if (f >= 0) {
g = (-1) * abs(sqrt(s)) -abs(sqrt(s))
} else { } else {
g = abs(sqrt(s)) abs(sqrt(s))
} }
val h = f * g - s val h = f * g - s
this[i, l] = f - g this[i, l] = f - g
@ -457,7 +462,7 @@ internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2
for (i in min(n, m) - 1 downTo 0) { for (i in min(n, m) - 1 downTo 0) {
l = i + 1 l = i + 1
g = wBuffer[wStart + i] g = wBuffer[wStart + i]
for (j in l until n) { for (j in l until n) {
this[i, j] = 0.0 this[i, j] = 0.0
} }
@ -484,13 +489,13 @@ internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2
this[i, i] += 1.0 this[i, i] += 1.0
} }
var flag = 0 var flag: Int
var nm = 0 var nm = 0
var c = 0.0 var c: Double
var h = 0.0 var h: Double
var y = 0.0 var y: Double
var z = 0.0 var z: Double
var x = 0.0 var x: Double
for (k in n - 1 downTo 0) { for (k in n - 1 downTo 0) {
for (its in 1 until iterations) { for (its in 1 until iterations) {
flag = 1 flag = 1
@ -531,7 +536,7 @@ internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2
} }
} }
z = wBuffer[wStart + k] z = wBuffer[wStart + k]
if (l == k) { if (l == k) {
if (z < 0.0) { if (z < 0.0) {
wBuffer[wStart + k] = -z wBuffer[wStart + k] = -z
@ -541,9 +546,9 @@ internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2
break break
} }
x = wBuffer[wStart + l] x = wBuffer[wStart + l]
nm = k - 1 nm = k - 1
y = wBuffer[wStart + nm] y = wBuffer[wStart + nm]
g = rv1[nm] g = rv1[nm]
h = rv1[k] h = rv1[k]
f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y) f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y)
@ -552,7 +557,7 @@ internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2
c = 1.0 c = 1.0
s = 1.0 s = 1.0
var i = 0 var i: Int
for (j in l until nm + 1) { for (j in l until nm + 1) {
i = j + 1 i = j + 1
g = rv1[i] g = rv1[i]