From 0b2b1b76e69e0c4edca5cf9209da42834b93d24f Mon Sep 17 00:00:00 2001
From: daniel <daniel@kapla.at>
Date: Mon, 16 Sep 2019 11:57:10 +0200
Subject: [PATCH] Improved runtime of pure R grad.

---
 CVE_R/R/gradient.R | 24 ++++++++++++------------
 test.R             |  2 +-
 wip.R              |  4 ++--
 3 files changed, 15 insertions(+), 15 deletions(-)

diff --git a/CVE_R/R/gradient.R b/CVE_R/R/gradient.R
index ed61645..f5cb257 100644
--- a/CVE_R/R/gradient.R
+++ b/CVE_R/R/gradient.R
@@ -40,15 +40,15 @@ grad <- function(X, Y, V, h,
     # Vectorized distance matrix `D`.
     vecD <- colSums(tcrossprod(Q, X_diff)^2)
 
-    # Weight matrix `W` (dnorm ... gaussean density function)
-    W <- matrix(1, n, n) # `exp(0) == 1`
-    W[lower] <- exp((-0.5 / h) * vecD^2) # Set lower tri. part
-    W[upper] <- t.default(W)[upper] # Mirror lower tri. to upper
-    W <- sweep(W, 2, colSums(W), FUN = `/`) # Col-Normalize
+    # Create Kernel matrix (aka. apply kernel to distances)
+    K <- matrix(1, n, n) # `exp(0) == 1`
+    K[lower] <- exp((-0.5 / h) * vecD^2) # Set lower tri. part
+    K[upper] <- t(K)[upper] # Mirror lower tri. to upper
 
     # Weighted `Y` momentums
-    y1 <- Y %*% W # Result is 1D -> transposition irrelevant
-    y2 <- Y^2 %*% W
+    colSumsK <- colSums(K)
+    y1 <- (K %*% Y) / colSumsK
+    y2 <- (K %*% Y^2) / colSumsK
 
     # Per example loss `L(V, X_i)`
     L <- y2 - y1^2
@@ -59,11 +59,11 @@ grad <- function(X, Y, V, h,
         loss <<- mean(L)
     }
 
-    # Vectorized Weights with forced symmetry
-    vecS <- (L[i] - (Y[j] - y1[i])^2) * W[lower]
-    vecS <- vecS + ((L[j] - (Y[i] - y1[j])^2) * W[upper])
-    # Compute scaling of `X` row differences
-    vecS <- vecS * vecD
+    # Compute scaling vector `vecS` for `X_diff`.
+    tmp <- kronecker(matrix(y1, n, 1), matrix(Y, 1, n), `-`)^2
+    tmp <- as.vector(L) - tmp
+    tmp <- tmp * K / colSumsK
+    vecS <- (tmp + t(tmp))[lower] * vecD
 
     # The gradient.
     # 1. The `crossprod(A, B)` is equivalent to `t(A) %*% B`,
diff --git a/test.R b/test.R
index 5fa890b..5826130 100644
--- a/test.R
+++ b/test.R
@@ -2,7 +2,7 @@
 # path <- '~/Projects/CVE/tmp/logger.R.pdf'
 
 library(CVE)
-path <- '~/Projects/CVE/tmp/seeded_test.pdf'
+path <- '~/Projects/CVE/tmp/logger.C.pdf'
 
 epochs <- 100
 attempts <- 25
diff --git a/wip.R b/wip.R
index 7648c37..03252b5 100644
--- a/wip.R
+++ b/wip.R
@@ -55,11 +55,10 @@ grad2 <- function(X, Y, V, h, persistent = TRUE) {
     # vecD <- rowSums((X_diff %*% Q)^2)
     vecD <- colSums(tcrossprod(Q, X_diff)^2)
 
-    # Weight matrix `W` (dnorm ... gaussean density function)
+    # Create Kernel matrix (aka. apply kernel to distances)
     K <- matrix(1, n, n) # `exp(0) == 1`
     K[lower] <- exp((-0.5 / h) * vecD^2) # Set lower tri. part
     K[upper] <- t(K)[upper] # Mirror lower tri. to upper
-    # W <- sweep(K, 2, colSums(K), FUN = `/`) # Col-Normalize
 
     # Weighted `Y` momentums
     colSumsK <- colSums(K)
@@ -68,6 +67,7 @@ grad2 <- function(X, Y, V, h, persistent = TRUE) {
     # Per example loss `L(V, X_i)`
     L <- y2 - y1^2
 
+    # Compute scaling vector `vecS` for `X_diff`.
     tmp <- kronecker(matrix(y1, n, 1), matrix(Y, 1, n), `-`)^2
     tmp <- as.vector(L) - tmp
     tmp <- tmp * K / colSumsK