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

Feedback for SVD

Browse Source
This commit is contained in:
Roland Grinis 2021-04-26 16:24:26 +01:00
parent 30ca333c04
commit c2db3a23e1
1 changed files with 10 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
examples/src/main/kotlin/space/kscience/kmath/tensors/OLSWithSVD.kt
View File

@ -9,6 +9,8 @@ import space.kscience.kmath.tensors.core.DoubleTensor
import space.kscience.kmath.tensors.core.algebras.DoubleAnalyticTensorAlgebra import space.kscience.kmath.tensors.core.algebras.DoubleAnalyticTensorAlgebra
import space.kscience.kmath.tensors.core.algebras.DoubleLinearOpsTensorAlgebra import space.kscience.kmath.tensors.core.algebras.DoubleLinearOpsTensorAlgebra
import kotlin.math.abs
// OLS estimator using SVD // OLS estimator using SVD
fun main() { fun main() {
@ -26,8 +28,7 @@ fun main() {
doubleArrayOf(1.0, 2.5, 3.4, 5.0, 10.1) doubleArrayOf(1.0, 2.5, 3.4, 5.0, 10.1)
) )
println("Real alpha:\n" + println("Real alpha:\n$alpha")
"$alpha")
// also take sample of size 20 from normal distribution for x TODO rename // also take sample of size 20 from normal distribution for x TODO rename
val x = randNormal( val x = randNormal(
@ -35,20 +36,22 @@ fun main() {
randSeed randSeed
) )
// calculate y and add gaussian noise (N(0, 0.05)) TODO rename // calculate y and add gaussian noise (N(0, 0.05))
// TODO: please add an intercept: Y = beta * X + alpha + N(0,0.5)
val y = x dot alpha val y = x dot alpha
y += y.randNormalLike(randSeed) * 0.05 y += y.randNormalLike(randSeed) * 0.05
// now restore the coefficient vector with OSL estimator with SVD // now restore the coefficient vector with OSL estimator with SVD
// TODO: you need to change accordingly [X 1] [alpha beta] = Y
// TODO: inverting [X 1] via SVD
val (u, singValues, v) = x.svd() val (u, singValues, v) = x.svd()
// we have to make sure the singular values of the matrix are not close to zero // we have to make sure the singular values of the matrix are not close to zero
println("Singular values:\n" + println("Singular values:\n$singValues")
"$singValues")
// TODO something with Boolean tensors
// inverse Sigma matrix can be restored from singular values with diagonalEmbedding function // inverse Sigma matrix can be restored from singular values with diagonalEmbedding function
val sigma = diagonalEmbedding(1.0/singValues) val sigma = diagonalEmbedding(singValues.map{ x -> if (abs(x) < 1e-3) 0.0 else 1.0/x })
val alphaOLS = v dot sigma dot u.transpose() dot y val alphaOLS = v dot sigma dot u.transpose() dot y
println("Estimated alpha:\n" + println("Estimated alpha:\n" +