work on regularization term with nonlinearity

3 files changed, 27 insertions, 15 deletions

diff --git a/src/msvmmaj_kernel.c b/src/msvmmaj_kernel.c
index fc699dd..5ac138c 100644
--- a/src/msvmmaj_kernel.c
+++ b/src/msvmmaj_kernel.c
@@ -31,14 +31,18 @@
  */
 void msvmmaj_make_kernel(struct MajModel *model, struct MajData *data)
 {
-	if (model->kerneltype == K_LINEAR)
+	long i, j;
+	if (model->kerneltype == K_LINEAR) {
+		model->J = Calloc(double, model->m+1);
+		for (i=1; i<model->m+1; i++) 
+			matrix_set(model->J, 1, i, 0, 1.0);
 		return;
+	}
 
-	long i, j;
 	long n = model->n;
 	double value;
 	double *x1, *x2;
-	double *K = Calloc(double, n*n*sizeof(double));
+	double *K = Calloc(double, n*n);
 
 	for (i=0; i<n; i++) {
 		for (j=i; j<n; j++) {
@@ -63,8 +67,6 @@ void msvmmaj_make_kernel(struct MajModel *model, struct MajData *data)
 		}
 	}
 
-	print_matrix(K, n, n);
-
 	double *P = Malloc(double, n*n);
 	double *Lambda = Malloc(double, n);
 	long num_eigen = msvmmaj_make_eigen(K, n, P, Lambda);
@@ -79,6 +81,13 @@ void msvmmaj_make_kernel(struct MajModel *model, struct MajData *data)
 	}
 	data->m = n;
 
+	// Set the regularization matrix (change if not full rank used)
+	model->J = Calloc(double, model->m+1);
+	for (i=1; i<model->m+1; i++) {
+		value = 1.0/matrix_get(Lambda, 1, i-1, 0);
+		matrix_set(model->J, 1, i, 0, value);
+	}
+
 	// let data know what it's made of
 	data->kerneltype = model->kerneltype;
 	free(data->kernelparam);
@@ -192,8 +201,6 @@ long msvmmaj_make_eigen(double *K, long n, double *P, double *Lambda)
 		for (j=0; j<n; j++)
 			P[i*n+j] = tempP[j*n+i];
 
-	print_matrix(P, n, n);
-
 	free(tempP);
 
 	// replace by number of columns of P
diff --git a/src/msvmmaj_train.c b/src/msvmmaj_train.c
index 97ee6a1..311b2d4 100644
--- a/src/msvmmaj_train.c
+++ b/src/msvmmaj_train.c
@@ -86,7 +86,7 @@ void msvmmaj_optimize(struct MajModel *model, struct MajData *data)
 		Lbar = L;
 		L = msvmmaj_get_loss(model, data, ZV);
 
-		if (it%50 == 0)
+		if (it%1 == 0)
 			note("iter = %li, L = %15.16f, Lbar = %15.16f, "
 			     "reldiff = %15.16f\n", it, L, Lbar, (Lbar - L)/L);
 		it++;
@@ -153,10 +153,12 @@ double msvmmaj_get_loss(struct MajModel *model, struct MajData *data,
 	loss /= ((double) n);
 
 	value = 0;
-	for (i=1; i<m+1; i++) {
+	for (i=0; i<m+1; i++) {
+		rowvalue = 0;
 		for (j=0; j<K-1; j++) {
-			value += pow(matrix_get(model->V, K-1, i, j), 2.0);
+			rowvalue += pow(matrix_get(model->V, K-1, i, j), 2.0);
 		}
+		value += model->J[i] * rowvalue;
 	}
 	loss += model->lambda * value;
 
@@ -418,12 +420,14 @@ void msvmmaj_get_update(struct MajModel *model, struct MajData *data, double *B,
 			ZAZV,
 			K-1);
 	/* 
-	 * Add lambda to all diagonal elements except the 
-	 * first one. 
+	 * Add lambda to all diagonal elements except the first one. Recall 
+	 * that ZAZ is of size m+1 and is symmetric.
 	 */
 	i = 0;
-	for (j=0; j<m; j++)
-		ZAZ[i+=m+1 + 1] += model->lambda;
+	for (j=0; j<m; j++) {
+		i += (m+1) + 1;
+		ZAZ[i] += model->lambda * model->J[j+1];
+	}
 	
 	// For the LAPACK call we need to switch to Column-
 	// Major order. This is unnecessary for the matrix
diff --git a/src/trainMSVMMaj.c b/src/trainMSVMMaj.c
index 66f6450..af91eaf 100644
--- a/src/trainMSVMMaj.c
+++ b/src/trainMSVMMaj.c
@@ -105,7 +105,8 @@ int main(int argc, char **argv)
 
 	// seed the random number generator (only place in programs is in
 	// command line interfaces)
-	srand(time(NULL));
+	//srand(time(NULL));
+	srand(123456);
 
 	if (msvmmaj_check_argv_eq(argc, argv, "-m")) {
 		struct MajModel *seed_model = msvmmaj_init_model();