Rewrite UU matrix to be K*K*(K-1) instead of n*K*(K-1)

significant memory reduction by turning the 3D UU matrix into a
2D block matrix, with significantly less dimensions

4 files changed, 55 insertions, 63 deletions

diff --git a/src/gensvm_base.c b/src/gensvm_base.c
index 2c30530..c77e704 100644
--- a/src/gensvm_base.c
+++ b/src/gensvm_base.c
@@ -119,7 +119,7 @@ void gensvm_allocate_model(struct GenModel *model)
 	model->V = Calloc(double, (m+1)*(K-1));
 	model->Vbar = Calloc(double, (m+1)*(K-1));
 	model->U = Calloc(double, K*(K-1));
-	model->UU = Calloc(double, n*K*(K-1));
+	model->UU = Calloc(double, K*K*(K-1));
 	model->Q = Calloc(double, n*K);
 	model->H = Calloc(double, n*K);
 	model->rho = Calloc(double, n);
@@ -145,9 +145,6 @@ void gensvm_reallocate_model(struct GenModel *model, long n, long m)
 	if (model->n == n && model->m == m)
 		return;
 	if (model->n != n) {
-		model->UU = Realloc(model->UU, double, n*K*(K-1));
-		Memset(model->UU, double, n*K*(K-1));
-
 		model->Q = Realloc(model->Q, double, n*K);
 		Memset(model->Q, double, n*K);
 
diff --git a/src/gensvm_optimize.c b/src/gensvm_optimize.c
index d08f917..a16512f 100644
--- a/src/gensvm_optimize.c
+++ b/src/gensvm_optimize.c
@@ -66,9 +66,8 @@ void gensvm_optimize(struct GenModel *model, struct GenData *data)
 	note("\tepsilon = %g\n", model->epsilon);
 	note("\n");
 
-	gensvm_simplex(model->K, model->U);
-	gensvm_simplex_diff(model, data);
 	gensvm_simplex(model);
+	gensvm_simplex_diff(model);
 
 	L = gensvm_get_loss(model, data, ZV);
 	Lbar = L + 2.0*model->epsilon*L;
@@ -576,39 +575,6 @@ void gensvm_get_update(struct GenModel *model, struct GenData *data,
 	}
 }
 
-/**
- * @brief Generate the simplex difference matrix
- *
- * @details
- * The simplex difference matrix is a 3D matrix which is constructed
- * as follows. For each instance i, the difference vectors between the row of
- * the simplex matrix corresponding to the class label of instance i and the
- * other rows of the simplex matrix are calculated. These difference vectors
- * are stored in a matrix, which is one horizontal slice of the 3D matrix.
- *
- * We use the indices i, j, k for the three dimensions n, K-1, K of UU. Then
- * the i,j,k -th element of UU is equal to U(y[i]-1, j) - U(k, j).
- *
- * @param[in,out] 	model 	the corresponding GenModel
- * @param[in] 		data 	the corresponding GenData
- *
- */
-void gensvm_simplex_diff(struct GenModel *model, struct GenData *data)
-{
-	long i, j, k;
-	double value;
-
-	long n = model->n;
-	long K = model->K;
-
-	for (i=0; i<n; i++) {
-		for (j=0; j<K-1; j++) {
-			for (k=0; k<K; k++) {
-				value = matrix_get(model->U, K-1,
-						data->y[i]-1, j);
-				value -= matrix_get(model->U, K-1, k, j);
-				matrix3_set(model->UU, K-1, K, i, j, k, value);
-			}
 		}
 	}
 }
@@ -689,21 +655,17 @@ void gensvm_calculate_huber(struct GenModel *model)
  * allocated. In addition, the matrix ZV is calculated here. It is assigned
  * to a pre-allocated block of memory, which is passed to this function.
  *
- * @todo
- * Transform UU to small UU then fix that here
- *
  * @param[in,out] 	model 	the corresponding GenModel
  * @param[in] 		data 	the corresponding GenData
  * @param[in,out] 	ZV 	a pointer to a memory block for ZV. On exit
  * 				this block is updated with the new ZV matrix
  * 				calculated with GenModel::V
- *
  */
 void gensvm_calculate_errors(struct GenModel *model, struct GenData *data,
 		double *ZV)
 {
-	long i, j, k;
-	double zv, value;
+	long i, j;
+	double q, *uu_row;
 
 	long n = model->n;
 	long m = model->m;
@@ -725,15 +687,13 @@ void gensvm_calculate_errors(struct GenModel *model, struct GenData *data,
 			ZV,
 			K-1);
 
-	Memset(model->Q, double, n*K);
 	for (i=0; i<n; i++) {
-		for (j=0; j<K-1; j++) {
-			zv = matrix_get(ZV, K-1, i, j);
-			for (k=0; k<K; k++) {
-				value = zv * matrix3_get(model->UU, K-1, K, i,
-						j, k);
-				matrix_add(model->Q, K, i, k, value);
-			}
+		for (j=0; j<K; j++) {
+			if (j == (data->y[i]-1))
+				continue;
+			uu_row = &model->UU[((data->y[i]-1)*K+j)*(K-1)];
+			q = cblas_ddot(K-1, &ZV[i*(K-1)], 1, uu_row, 1);
+			matrix_set(model->Q, K, i, j, q);
 		}
 	}
 }
diff --git a/src/gensvm_pred.c b/src/gensvm_pred.c
index 8a9a43e..43b27cc 100644
--- a/src/gensvm_pred.c
+++ b/src/gensvm_pred.c
@@ -31,7 +31,7 @@ void gensvm_predict_labels(struct GenData *testdata, struct GenModel *model,
 		long *predy)
 {
 	long i, j, k, n, m, K, label;
-	double norm, min_dist, *S, *ZV, *U;
+	double norm, min_dist, *S, *ZV;
 
 	n = testdata->n;
 	m = testdata->r;
@@ -40,10 +40,9 @@ void gensvm_predict_labels(struct GenData *testdata, struct GenModel *model,
 	// allocate necessary memory
 	S = Calloc(double, K-1);
 	ZV = Calloc(double, n*(K-1));
-	U = Calloc(double, K*(K-1));
 
 	// Generate the simplex matrix
-	gensvm_simplex(K, U);
+	gensvm_simplex(model);
 
 	// Generate the simplex space vectors
 	cblas_dgemm(
@@ -70,7 +69,7 @@ void gensvm_predict_labels(struct GenData *testdata, struct GenModel *model,
 		for (j=0; j<K; j++) {
 			for (k=0; k<K-1; k++) {
 				S[k] = matrix_get(ZV, K-1, i, k) -
-					matrix_get(U, K-1, j, k);
+					matrix_get(model->U, K-1, j, k);
 			}
 			norm = cblas_dnrm2(K-1, S, 1);
 			if (norm < min_dist) {
@@ -82,7 +81,6 @@ void gensvm_predict_labels(struct GenData *testdata, struct GenModel *model,
 	}
 
 	free(ZV);
-	free(U);
 	free(S);
 }
 
diff --git a/src/gensvm_simplex.c b/src/gensvm_simplex.c
index 1fd5f14..a704d85 100644
--- a/src/gensvm_simplex.c
+++ b/src/gensvm_simplex.c
@@ -25,21 +25,58 @@
  * @param[in] 		K 	number of classes
  * @param[in,out] 	U 	simplex matrix of size K * (K-1)
  */
-void gensvm_simplex(long K, double *U)
+void gensvm_simplex(struct GenModel *model)
 {
-	long i, j;
+	long i, j, K = model->K;
+
 	for (i=0; i<K; i++) {
 		for (j=0; j<K-1; j++) {
 			if (i <= j) {
-				matrix_set(U, K-1, i, j,
+				matrix_set(model->U, K-1, i, j,
 					       	-1.0/sqrt(2.0*(j+1)*(j+2)));
 			} else if (i == j+1) {
-				matrix_set(U, K-1, i, j,
+				matrix_set(model->U, K-1, i, j,
 					       	sqrt((j+1)/(2.0*(j+2))));
 			} else {
-				matrix_set(U, K-1, i, j, 0.0);
+				matrix_set(model->U, K-1, i, j, 0.0);
 			}
 		}
 	}
 }
 
+/**
+ * @brief Generate the simplex difference matrix
+ *
+ * @details
+ * The simplex difference matrix is a 2D block matrix which is constructed
+ * as follows. For each class i, we have a block of K rows and K-1 columns.
+ * Each row in the block for class i contains a row vector with the difference
+ * of the simplex matrix, U(i, :) - U(j, :).
+ *
+ * In the paper the notation @f$\boldsymbol{\delta}_{kj}'@f$ is used for the
+ * difference vector of @f$\textbf{u}_k' - \textbf{u}_j'@f$, where
+ * @f$\textbf{u}_k'@f$ corresponds to row k of @f$\textbf{U}@f$. Due to the
+ * indexing in the paper being 1-based and C indexing is 0 based, the vector
+ * @f$\boldsymbol{\delta}_{kj}'@f$ corresponds to the row (k-1)*K+(j-1) in the
+ * UU matrix generated here.
+ *
+ * @param[in,out] 	model 	the corresponding GenModel
+ *
+ */
+void gensvm_simplex_diff(struct GenModel *model)
+{
+	long i, j, l, K = model->K;
+	double value;
+
+	// UU is a 2D block matrix, where block i has the differences:
+	// U(i, :) - U(j, :) for all j
+	for (i=0; i<K; i++) {
+		for (j=0; j<K; j++) {
+			for (l=0; l<K-1; l++) {
+				value = matrix_get(model->U, K-1, i, l);
+				value -= matrix_get(model->U, K-1, j, l);
+				matrix_set(model->UU, K-1, i*K+j, l, value);
+			}
+		}
+	}
+}