fiddling with cholesky backprop

GPy/util/choleskies.py

@@ -100,7 +100,7 @@ def indexes_to_fix_for_low_rank(rank, size):
 if config.getboolean('cython', 'working'):
     triang_to_flat = _triang_to_flat_cython
     flat_to_triang = _flat_to_triang_cython
-    backprop_gradient = choleskies_cython.backprop_gradient
+    backprop_gradient = choleskies_cython.backprop_gradient_par_c
 else:
     backprop_gradient = _backprop_gradient_pure
     triang_to_flat =  _triang_to_flat_pure

GPy/util/choleskies_cython.pyx

@@ -76,12 +76,36 @@ def backprop_gradient_par(double[:,:] dL, double[:,:] L):
         dL_dK[k, k] /= (2. * L[k, k])
     return dL_dK
 
-cdef extern from "cholesky_backprop.h":
+#here's a pure C version...
+cdef extern from "cholesky_backprop.h" nogil:
     void chol_backprop(int N, double* dL, double* L)
 
 def backprop_gradient_par_c(np.ndarray[double, ndim=2] dL, np.ndarray[double, ndim=2] L):
     cdef np.ndarray[double, ndim=2] dL_dK = np.tril(dL) # makes a copy, c-contig
     cdef int N = L.shape[0]
-    chol_backprop(N, <double*> dL_dK.data, <double*> L.data)
+    with nogil:
+        chol_backprop(N, <double*> dL_dK.data, <double*> L.data)
     return dL_dK
 
+
+# TODO: with the next release of cython, cimport scipy.linalg.cython_blas as blas, then blas the hell out of this.
+def backprop_gradient_par2(double[:,:] dL, double[:,:] L):
+    """
+    a very slow implementation, but the clearest I hope
+    """
+    cdef double[:,:] dL_dK = np.tril(dL).copy()
+    cdef int N = L.shape[0]
+    cdef int k, j, i, iN, kN
+    for k in range(N - 1, -1, -1):
+        #pragma this loop:
+        for i in range(k+1, N):
+            dL_dK[i,+k] -= np.dot(dL_dK[i,k+1:], L[k+1:,k])
+            dL_dK[i,+k] -= np.dot(dL_dK[i:,i], L[i:,k])
+
+        for i in range(k+1, N):
+            dL_dK[i, k] /= L[k, k];
+            dL_dK[k, k] -= L[i, k] * dL_dK[i, k];
+
+        dL_dK[k, k] /= (2. * L[k, k]);
+    return np.asarray(dL_dK)
+

GPy/util/cholesky_backprop.c

@@ -1,22 +1,32 @@
+#include <omp.h>
+double mydot(int n, double* a,  int stride_a, double* b, int stride_b){
+    double ret = 0;	
+    for(int i=0; i<n; i++){
+	ret += a[i*stride_a]*b[i*stride_b];
+    }
+    return ret;
+}
 
-void chol_backprop(int N, double* dL, double* L){
+
+void chol_backprop(int N, double* dL, double* U){
     //at the input to this fn, dL is df_dL. after this fn is complet, dL is df_dK
-    int i,j,k;
-    for(k=N-1;k>(-1);k--){
-        #pragma omp parallel for private(i,j)
-        for(i=k+1;i<N; i++){
-            for(j=k+1;j<(i+1);j++){
-                dL[i*N + k] -= dL[i *N + j] * L[j*N + k];
-	    }
-            for(j=i;j<N;j++){
-                dL[i*N + k] -= dL[j*N + i] * L[j*N +k];
-	    }
+    int iN, kN;
+    for(int k=N-1;k>(-1);k--){
+	kN = k*N;
+        #pragma omp parallel for private(iN)
+        for(int i=k+1; i<N; i++){
+	    iN = i*N;
+	    dL[iN+k] -= mydot(i-k, &dL[iN+k+1], 1, &U[kN+k+1], 1);
+	    dL[iN+k] -= mydot(N-i, &dL[iN+i], N, &U[kN+i],  1);
+
 	}
-        for(i=k + 1; i<N; i++){
-            dL[i*N + k] /= L[k*N + k];
-            dL[k*N + k] -= L[i*N + k] * dL[i*N + k];
+        for(int i=(k + 1); i<N; i++){
+	    iN = i*N;
+            dL[iN + k] /= U[kN + k];
+            dL[kN + k] -= U[kN + i] * dL[iN + k];
 	}
-        dL[k*N + k] /= (2. * L[k*N + k]);
+
+        dL[kN + k] /= (2. * U[kN + k]);
     }
 }
 

GPy/util/cholesky_backprop.h

@@ -1 +1,2 @@
+double mydot(int n, double* a,  int stride_a, double* b, int stride_b);
 void chol_backprop(int N, double* dL, double* L);