Parameter D in Eigenbasis Calculation from Pivoted Cholesky Output

src/main/scala/scalismo/kernels/Kernel.scala

@@ -273,7 +273,7 @@ object Kernel {
 
     // we compute the eigenvectors only approximately, to a tolerance of 1e-5. As the nystrom approximation is
     // anyway not exact, this should be sufficient for all practical cases.
-    val (uMat, lambdaMat) = PivotedCholesky.computeApproximateEig(k, ptsForNystrom, 1.0, RelativeTolerance(1e-5))
+    val (uMat, lambdaMat) = PivotedCholesky.computeApproximateEig(k, ptsForNystrom, RelativeTolerance(1e-5))
 
     val lambda = lambdaMat.map(lmbda => (lmbda / effectiveNumberOfPointsSampled.toDouble))
     val numParams = (for (i <- (0 until lambda.size) if lambda(i) >= 1e-8) yield 1).size

src/main/scala/scalismo/numerics/PivotedCholesky.scala

@@ -183,12 +183,12 @@ object PivotedCholesky {
     computeApproximateCholeskyGeneric[Int](kernel, indices, stoppingCriterion)
   }
 
-  private def computeApproximateEigGeneric[A](k: (A, A) => Double, xs: IndexedSeq[A], D: Double, sc: StoppingCriterion) = {
+  private def computeApproximateEigGeneric[A](k: (A, A) => Double, xs: IndexedSeq[A], sc: StoppingCriterion) = {
 
     val PivotedCholesky(l, _, _) = computeApproximateCholeskyGeneric(k, xs, sc)
 
-    val LD = l(::, 0 until l.cols).t *:* D
-    val phi: DenseMatrix[Double] = LD * l(::, 0 until l.cols)
+    val Lt = l(::, 0 until l.cols).t
+    val phi: DenseMatrix[Double] = Lt * l(::, 0 until l.cols)
 
     val SVD(v, _, _) = breeze.linalg.svd(phi)
     val U: DenseMatrix[Double] = l(::, 0 until l.cols) * v
@@ -209,23 +209,22 @@ object PivotedCholesky {
 
   }
 
-  def computeApproximateEig(A: DenseMatrix[Double], D: Double, sc: StoppingCriterion) = {
+  def computeApproximateEig(A: DenseMatrix[Double], sc: StoppingCriterion) = {
     val kernel: (Int, Int) => Double = (i, j) => A(i, j)
     val indices = IndexedSeq.range(0, A.cols)
 
-    extractEigenvalues(computeApproximateEigGeneric(kernel, indices, D, sc))
+    extractEigenvalues(computeApproximateEigGeneric(kernel, indices, sc))
   }
 
   def computeApproximateEig[D: NDSpace](kernel: MatrixValuedPDKernel[D],
-    xs: IndexedSeq[Point[D]], D: Double,
-    stoppingCriterion: StoppingCriterion) = {
+    xs: IndexedSeq[Point[D]], stoppingCriterion: StoppingCriterion) = {
 
     case class PointWithDim(point: Point[D], dim: Int)
     val dim = kernel.outputDim
     val xsWithDim: IndexedSeq[PointWithDim] = xs.flatMap(f => (0 until dim).map(i => PointWithDim(f, i)))
     def kscalar(x: PointWithDim, y: PointWithDim): Double = kernel(x.point, y.point)(x.dim, y.dim)
 
-    extractEigenvalues(computeApproximateEigGeneric[PointWithDim](kscalar, xsWithDim, D, stoppingCriterion))
+    extractEigenvalues(computeApproximateEigGeneric[PointWithDim](kscalar, xsWithDim, stoppingCriterion))
 
   }
 
@@ -236,7 +235,7 @@ object PivotedCholesky {
 
     def k(x: Point[D], y: Point[D]): Double = kernel(x, y)
 
-    extractEigenvalues(computeApproximateEigGeneric[Point[D]](k, xs, D, stoppingCriterion))
+    extractEigenvalues(computeApproximateEigGeneric[Point[D]](k, xs, stoppingCriterion))
 
   }
 }

src/main/scala/scalismo/statisticalmodel/DiscreteGaussianProcess.scala

@@ -122,7 +122,6 @@ class DiscreteGaussianProcess[D: NDSpace, +DDomain <: DiscreteDomain[D], Value]
 
     val (basis, scale) = PivotedCholesky.computeApproximateEig(
       cov.asBreezeMatrix,
-      D = 1.0,
       RelativeTolerance(0.0)
     )
 

src/main/scala/scalismo/statisticalmodel/DiscreteLowRankGaussianProcess.scala

@@ -450,11 +450,11 @@ object DiscreteLowRankGaussianProcess {
     // Whether it is more efficient to compute the PCA using the gram matrix or the covariance matrix, depends on
     // whether we have more examples than variables.
     val (basisMat, varianceVector) = if (X.rows > X.cols) {
-      val (u, d2) = PivotedCholesky.computeApproximateEig(X.t * X * (1.0 / (n - 1)), 1.0, stoppingCriterion)
+      val (u, d2) = PivotedCholesky.computeApproximateEig(X.t * X * (1.0 / (n - 1)), stoppingCriterion)
       (u, d2)
     } else {
 
-      val (v, d2) = PivotedCholesky.computeApproximateEig(X * X.t * (1.0 / (n - 1)), 1.0, stoppingCriterion)
+      val (v, d2) = PivotedCholesky.computeApproximateEig(X * X.t * (1.0 / (n - 1)), stoppingCriterion)
 
       val D = d2.map(v => Math.sqrt(v))
       val Dinv = D.map(d => if (d > 1e-6) 1.0 / d else 0.0)

src/main/scala/scalismo/statisticalmodel/LowRankGaussianProcess.scala

@@ -266,7 +266,6 @@ object LowRankGaussianProcess {
     val (basis, scale) = PivotedCholesky.computeApproximateEig(
       kernel = gp.cov,
       xs = domain.points.toIndexedSeq,
-      D = 1.0,
       stoppingCriterion = RelativeTolerance(relativeTolerance)
     )
 

src/test/scala/scalismo/numerics/PivotedCholeskyTest.scala

@@ -35,7 +35,7 @@ class PivotedCholeskyTest extends ScalismoTestSuite {
       val k = GaussianKernel[_1D](1.0)
       val matrixValuedK = DiagonalKernel[_1D](k, 1)
       val m = Kernel.computeKernelMatrix[_1D](pts, matrixValuedK)
-      val eigCholesky = PivotedCholesky.computeApproximateEig(matrixValuedK, pts, 1.0, PivotedCholesky.RelativeTolerance(1e-15))
+      val eigCholesky = PivotedCholesky.computeApproximateEig(matrixValuedK, pts, PivotedCholesky.RelativeTolerance(1e-15))
       val (u, d) = eigCholesky
       val D = (u * breeze.linalg.diag(d) * u.t) - m
       Math.sqrt(breeze.linalg.trace(D * D.t)) should be <= 1e-5
@@ -49,7 +49,7 @@ class PivotedCholeskyTest extends ScalismoTestSuite {
       val k = GaussianKernel[_3D](1.0)
       val matrixValuedK = DiagonalKernel[_3D](k, 3)
       val m = Kernel.computeKernelMatrix[_3D](pts, matrixValuedK)
-      val eigCholesky = PivotedCholesky.computeApproximateEig(matrixValuedK, pts, 1.0, PivotedCholesky.RelativeTolerance(1e-15))
+      val eigCholesky = PivotedCholesky.computeApproximateEig(matrixValuedK, pts, PivotedCholesky.RelativeTolerance(1e-15))
       val (u, d) = eigCholesky
       val D = (u * breeze.linalg.diag(d) * u.t) - m
       Math.sqrt(breeze.linalg.trace(D * D.t)) should be <= 1e-5