From 002f359382eeccbf7752fd95e178979f3243eb30 Mon Sep 17 00:00:00 2001 From: tekhnofiend Date: Fri, 18 Jan 2008 04:12:35 +0000 Subject: [PATCH] updated niche/fastica --- fastlib/u/niche/fastica/fastica.cc | 175 ++++++--------------- fastlib/u/niche/fastica/lin_alg.h | 243 +++++++++++++++++++++++------ 2 files changed, 246 insertions(+), 172 deletions(-) diff --git a/fastlib/u/niche/fastica/fastica.cc b/fastlib/u/niche/fastica/fastica.cc index b54e561aaf..ccec926ce9 100644 --- a/fastlib/u/niche/fastica/fastica.cc +++ b/fastlib/u/niche/fastica/fastica.cc @@ -23,85 +23,8 @@ namespace { } - void RandVector(Vector &v) { - - index_t d = v.length(); - v.SetZero(); - - for(index_t i = 0; i+1 < d; i+=2) { - double a = drand48(); - double b = drand48(); - double first_term = sqrt(-2 * log(a)); - double second_term = 2 * M_PI * b; - v[i] = first_term * cos(second_term); - v[i+1] = first_term * sin(second_term); - } - - if((d % 2) == 1) { - v[d - 1] = sqrt(-2 * log(drand48())) * cos(2 * M_PI * drand48()); - } - - la::Scale(1/sqrt(la::Dot(v, v)), &v); - - } - void Center(Matrix X, Matrix &X_centered) { - Vector col_vector_sum; - col_vector_sum.Init(X.n_rows()); - col_vector_sum.SetZero(); - - index_t n = X.n_cols(); - - for(index_t i = 0; i < n; i++) { - Vector cur_col_vector; - X.MakeColumnVector(i, &cur_col_vector); - la::AddTo(cur_col_vector, &col_vector_sum); - } - - la::Scale(1/(double) n, &col_vector_sum); - - X_centered.CopyValues(X); - - for(index_t i = 0; i < n; i++) { - Vector cur_col_vector; - X_centered.MakeColumnVector(i, &cur_col_vector); - la::SubFrom(col_vector_sum, &cur_col_vector); - } - - } - - - void Whiten(Matrix X, Matrix &X_whitened, Matrix &whitening_matrix, Matrix &dewhitening_matrix) { - Matrix X_cov, D, D_inv, E; - Vector D_vector; - - Scale(1 / (double) (X.n_cols() - 1), - MulTransBInit(&X, &X, &X_cov)); - - - la::EigenvectorsInit(X_cov, &D_vector, &E); - - E.set(0, 1, -E.get(0, 1)); - E.set(1, 1, -E.get(1, 1)); - - - - index_t d = D_vector.length(); - D.Init(d, d); - D.SetZero(); - D_inv.Init(d, d); - D_inv.SetZero(); - for(index_t i = 0; i < d; i++) { - double sqrt_val = sqrt(D_vector[i]); - D.set(i, i, sqrt_val); - D_inv.set(i, i, 1 / sqrt_val); - } - - la::MulTransBInit(D_inv, E, &whitening_matrix); - la::MulInit(E, D, &dewhitening_matrix); - la::MulInit(whitening_matrix, X, &X_whitened); - } void Orthogonalize(Matrix W_old, Matrix &W) { @@ -157,11 +80,11 @@ namespace { int FixedPointICA(Matrix X, Matrix whitening_matrix, Matrix dewhitening_matrix, - int approach, index_t num_of_IC, int g, int fine_tune, - double a1, double a2, double mu, bool stabilization, - double epsilon, index_t max_num_iterations, - index_t max_fine_tune, double sample_size, - Matrix* A, Matrix* W) { + int approach, index_t num_of_IC, int g, int fine_tune, + double a1, double a2, double mu, bool stabilization, + double epsilon, index_t max_num_iterations, + index_t max_fine_tune, double sample_size, + Matrix* A, Matrix* W) { // ensure default values are passed into this function if the user doesn't care about certain parameters @@ -258,11 +181,11 @@ namespace { B.Init(d, num_of_IC); /* - for(index_t i = 0; i < num_of_IC; i++) { + for(index_t i = 0; i < num_of_IC; i++) { Vector b; B.MakeColumnVector(i, &b); RandVector(b); - } + } */ B.SetZero(); for(index_t i = 0; i < d; i++) { @@ -393,12 +316,12 @@ namespace { Scale(1 / (double) n, AddTo(MulInit(&X, &hyp_tan, &temp1), - DotMultiplyOverwrite(MulInit(&col_vector, - MapOverwrite(&MinusArg, - n, - MatrixMapSum(&Square, 0, &hyp_tan, &sum)), - &temp2), - &B))); + DotMultiplyOverwrite(MulInit(&col_vector, + MapOverwrite(&MinusArg, + n, + MatrixMapSum(&Square, 0, &hyp_tan, &sum)), + &temp2), + &B))); break; @@ -453,12 +376,12 @@ namespace { Scale(1 / (double) num_selected, AddTo(MulInit(&X_sub, &hyp_tan, &temp1), - DotMultiplyOverwrite(MulInit(&col_vector, - MapOverwrite(&MinusArg, - num_selected, - MatrixMapSum(&Square, 0, &hyp_tan, &sum)), - &temp2), - &B))); + DotMultiplyOverwrite(MulInit(&col_vector, + MapOverwrite(&MinusArg, + num_selected, + MatrixMapSum(&Square, 0, &hyp_tan, &sum)), + &temp2), + &B))); break; @@ -607,11 +530,11 @@ namespace { Scale(1 / (double) num_selected, SubOverwrite(MulInit(&X_sub, &U, &temp1), - DotMultiplyOverwrite(&B, - MulInit(&col_vector, - Sum(&ex, &sum), - &temp2)), - &B)); + DotMultiplyOverwrite(&B, + MulInit(&col_vector, + Sum(&ex, &sum), + &temp2)), + &B)); break; @@ -660,8 +583,8 @@ namespace { //D = diag(1 ./ (Beta + sum((Y_squared_a2 - 1) .* ex))) VectorToDiag(MapOverwrite(&Inv, - 0, - AddTo(&Beta_vector, &sum_vector)), + 0, + AddTo(&Beta_vector, &sum_vector)), &D); //B = B + myy * B * (Y' * gauss - diag(Beta)) * D; @@ -734,8 +657,8 @@ namespace { Scale(1 / (double) num_selected, MulInit(&X_sub, MapOverwrite(&pow, - 3, - MulTransAInit(&X_sub, &B, &temp1)), + 3, + MulTransAInit(&X_sub, &B, &temp1)), &temp2)); AddTo(&temp2, @@ -784,10 +707,10 @@ namespace { Scale(1 / (double) n, MulOverwrite(&X, - MapOverwrite(&Square, - 0, - MulTransAInit(&X, &B, &temp1)), - &B)); + MapOverwrite(&Square, + 0, + MulTransAInit(&X, &B, &temp1)), + &B)); break; } @@ -824,10 +747,10 @@ namespace { Scale(1 / (double) num_selected, MulOverwrite(&X_sub, - MapOverwrite(&Square, - 0, - MulTransAInit(&X_sub, &B, &temp1)), - &B)); + MapOverwrite(&Square, + 0, + MulTransAInit(&X_sub, &B, &temp1)), + &B)); break; } @@ -1405,11 +1328,11 @@ namespace { } -int FastICAMain(datanode *module) { +int FastICAMain(datanode *module, Matrix *W, Matrix *Y) { const char *data = fx_param_str_req(NULL, "data"); - Matrix X, X_centered, X_whitened, whitening_matrix, dewhitening_matrix, A, W; + Matrix X, X_centered, X_whitened, whitening_matrix, dewhitening_matrix, A; data::Load(data, &X); index_t d = X.n_rows(); // number of dimensions @@ -1490,27 +1413,22 @@ int FastICAMain(datanode *module) { fx_timer_start(NULL, "FastICA"); - X_centered.Init(d, n); - Center(X, X_centered); + Center(X, &X_centered); - Whiten(X_centered, X_whitened, whitening_matrix, dewhitening_matrix); + WhitenUsingEig(X_centered, &X_whitened, &whitening_matrix, &dewhitening_matrix); - FixedPointICA(X_whitened, whitening_matrix, dewhitening_matrix, approach, num_of_IC, nonlinearity, fine_tune, a1, a2, mu, stabilization, epsilon, max_num_iterations, max_fine_tune, sample_size, &A, &W); + FixedPointICA(X_whitened, whitening_matrix, dewhitening_matrix, approach, num_of_IC, nonlinearity, fine_tune, a1, a2, mu, stabilization, epsilon, max_num_iterations, max_fine_tune, sample_size, &A, W); - W.PrintDebug("W"); + W -> PrintDebug("W"); - Matrix Y; - la::MulInit(W, X, &Y); + la::MulInit(*W, X, Y); fx_timer_stop(NULL, "FastICA"); - SaveCorrectly("unmixing_matrix.dat", W); - SaveCorrectly("indep_comps.dat", Y); - return SUCCESS_PASS; } @@ -1520,7 +1438,12 @@ int main(int argc, char *argv[]) { srand48(time(0)); - int ret_val = FastICAMain(fx_root); + Matrix W, Y; + + int ret_val = FastICAMain(fx_root, &W, &Y); + + SaveCorrectly("unmixing_matrix.dat", W); + SaveCorrectly("indep_comps.dat", Y); fx_done(); diff --git a/fastlib/u/niche/fastica/lin_alg.h b/fastlib/u/niche/fastica/lin_alg.h index ea9883cb43..d9f915e96d 100644 --- a/fastlib/u/niche/fastica/lin_alg.h +++ b/fastlib/u/niche/fastica/lin_alg.h @@ -37,6 +37,15 @@ double ArgMinus(double x, double arg) { return arg - x; } +Matrix* DiagMatrixInit(index_t n, double value, Matrix *diag_matrix) { + diag_matrix -> Init(n, n); + diag_matrix -> SetZero(); + for(index_t i = 0; i < n; i++) { + diag_matrix -> set(i, i, value); + } + + return diag_matrix; +} Matrix* ColVector(index_t n, double value, Matrix *col_vector) { col_vector -> Init(n, 1); @@ -490,68 +499,210 @@ void MakeSubMatrixByColumns(Vector column_indices, Matrix A, Matrix *A_sub) { } } + + +void Center(Matrix X, Matrix* X_centered) { + Vector col_vector_sum; + col_vector_sum.Init(X.n_rows()); + col_vector_sum.SetZero(); + + index_t n = X.n_cols(); + + for(index_t i = 0; i < n; i++) { + Vector cur_col_vector; + X.MakeColumnVector(i, &cur_col_vector); + la::AddTo(cur_col_vector, &col_vector_sum); + } + + la::Scale(1/(double) n, &col_vector_sum); + + X_centered -> Copy(X); + + for(index_t i = 0; i < n; i++) { + Vector cur_col_vector; + X_centered -> MakeColumnVector(i, &cur_col_vector); + la::SubFrom(col_vector_sum, &cur_col_vector); + } + +} + + +void WhitenUsingSVD(Matrix X, Matrix* X_whitened, Matrix* whitening_matrix) { + + Matrix cov_X, U, VT, inv_S_matrix, temp1; + Vector S_vector; + + Scale(1 / (double) (X.n_cols() - 1), + MulTransBInit(&X, &X, &cov_X)); + + la::SVDInit(cov_X, &S_vector, &U, &VT); + + index_t d = S_vector.length(); + inv_S_matrix.Init(d, d); + inv_S_matrix.SetZero(); + for(index_t i = 0; i < d; i++) { + double inv_sqrt_val = 1 / sqrt(S_vector[i]); + inv_S_matrix.set(i, i, inv_sqrt_val); + } + + cov_X.PrintDebug("cov(X')"); + U.PrintDebug("U"); + VT.PrintDebug("VT"); + inv_S_matrix.PrintDebug("S^-.5"); + + MulTransBInit(MulTransAInit(&VT, &inv_S_matrix, &temp1), + &U, + whitening_matrix); + + MulInit(whitening_matrix, &X, X_whitened); + +} + + +void WhitenUsingEig(Matrix X, Matrix* X_whitened, Matrix* whitening_matrix, Matrix* dewhitening_matrix) { + Matrix cov_X, D, D_inv, E; + Vector D_vector; + + Scale(1 / (double) (X.n_cols() - 1), + MulTransBInit(&X, &X, &cov_X)); + + + la::EigenvectorsInit(cov_X, &D_vector, &E); + + //E.set(0, 1, -E.get(0, 1)); + //E.set(1, 1, -E.get(1, 1)); + + + + index_t d = D_vector.length(); + D.Init(d, d); + D.SetZero(); + D_inv.Init(d, d); + D_inv.SetZero(); + for(index_t i = 0; i < d; i++) { + double sqrt_val = sqrt(D_vector[i]); + D.set(i, i, sqrt_val); + D_inv.set(i, i, 1 / sqrt_val); + } + + la::MulTransBInit(D_inv, E, whitening_matrix); + la::MulInit(E, D, dewhitening_matrix); + la::MulInit(*whitening_matrix, X, X_whitened); +} + + + +void RandVector(Vector &v) { + + index_t d = v.length(); + v.SetZero(); + + for(index_t i = 0; i+1 < d; i+=2) { + double a = drand48(); + double b = drand48(); + double first_term = sqrt(-2 * log(a)); + double second_term = 2 * M_PI * b; + v[i] = first_term * cos(second_term); + v[i+1] = first_term * sin(second_term); + } + + if((d % 2) == 1) { + v[d - 1] = sqrt(-2 * log(drand48())) * cos(2 * M_PI * drand48()); + } + + la::Scale(1/sqrt(la::Dot(v, v)), &v); + +} + + +Matrix* RandNormalInit(index_t d, index_t n, Matrix* A) { + + double* A_elements = A -> ptr(); + + index_t num_elements = d * n; + + for(index_t i = 0; i+1 < num_elements; i+=2) { + double a = drand48(); + double b = drand48(); + double first_term = sqrt(-2 * log(a)); + double second_term = 2 * M_PI * b; + A_elements[i] = first_term * cos(second_term); + A_elements[i+1] = first_term * sin(second_term); + } + + if((d % 2) == 1) { + A_elements[d - 1] = sqrt(-2 * log(drand48())) * cos(2 * M_PI * drand48()); + } + + return A; +} + + + + + /* int main(int argc, char *argv[]) {*/ - //fx_init(argc, argv); - /* - Matrix A, B, C, D, E, F; +//fx_init(argc, argv); +/* + Matrix A, B, C, D, E, F; - RandMatrix(5, 2, &A); - RandMatrix(2, 4, &B); - RandMatrix(5, 4, &C); - RandMatrix(7, 5, &D); - RandMatrix(7, 4, &E); - RandMatrix(7, 4, &F); + RandMatrix(5, 2, &A); + RandMatrix(2, 4, &B); + RandMatrix(5, 4, &C); + RandMatrix(7, 5, &D); + RandMatrix(7, 4, &E); + RandMatrix(7, 4, &F); - SaveCorrectly("A.dat", A); - SaveCorrectly("B.dat", B); - SaveCorrectly("C.dat", C); - SaveCorrectly("D.dat", D); - SaveCorrectly("E.dat", E); + SaveCorrectly("A.dat", A); + SaveCorrectly("B.dat", B); + SaveCorrectly("C.dat", C); + SaveCorrectly("D.dat", D); + SaveCorrectly("E.dat", E); - Matrix temp1, temp2, temp3, temp4, temp5; + Matrix temp1, temp2, temp3, temp4, temp5; - Vector sum_vector; - Matrix diag_matrix; + Vector sum_vector; + Matrix diag_matrix; - VectorToDiag(MatrixSum(DotMultiplyInit(MapOverwrite(&TimesTen, Sub(Mul(&D, Sub(Mul(&A, &B, &temp1), &C, &temp2), &temp3), &E, &temp4)), &F, &temp5), &sum_vector), &diag_matrix); + VectorToDiag(MatrixSum(DotMultiplyInit(MapOverwrite(&TimesTen, Sub(Mul(&D, Sub(Mul(&A, &B, &temp1), &C, &temp2), &temp3), &E, &temp4)), &F, &temp5), &sum_vector), &diag_matrix); - Vector diag_vector; - DiagToVector(&diag_matrix, &diag_vector); + Vector diag_vector; + DiagToVector(&diag_matrix, &diag_vector); - la::Scale(2, Scale(100, &diag_matrix)); + la::Scale(2, Scale(100, &diag_matrix)); - A.PrintDebug("A"); - B.PrintDebug("B"); - C.PrintDebug("C"); - D.PrintDebug("D"); - E.PrintDebug("E"); - F.PrintDebug("F"); + A.PrintDebug("A"); + B.PrintDebug("B"); + C.PrintDebug("C"); + D.PrintDebug("D"); + E.PrintDebug("E"); + F.PrintDebug("F"); - temp4.PrintDebug("temp4"); - temp5.PrintDebug("temp5"); + temp4.PrintDebug("temp4"); + temp5.PrintDebug("temp5"); - sum_vector.PrintDebug("sum(temp5)"); + sum_vector.PrintDebug("sum(temp5)"); - diag_matrix.PrintDebug("diag_matrix"); - diag_vector.PrintDebug("diag_vector"); + diag_matrix.PrintDebug("diag_matrix"); + diag_vector.PrintDebug("diag_vector"); - Matrix Z, sub_Z; - RandMatrix(4, 5, &Z); - Vector indices; - indices.Init(3); - indices[0] = 1; - indices[1] = 0; - indices[2] = 3; - MakeSubMatrixByColumns(indices, Z, &sub_Z); + Matrix Z, sub_Z; + RandMatrix(4, 5, &Z); + Vector indices; + indices.Init(3); + indices[0] = 1; + indices[1] = 0; + indices[2] = 3; + MakeSubMatrixByColumns(indices, Z, &sub_Z); - Z.PrintDebug("Z"); - sub_Z.PrintDebug("sub_Z"); + Z.PrintDebug("Z"); + sub_Z.PrintDebug("sub_Z"); - */ +*/ /* index_t n = 10; index_t d = 2; @@ -568,8 +719,8 @@ void MakeSubMatrixByColumns(Vector column_indices, Matrix A, Matrix *A_sub) { Matrix sum, temp1, temp2; AddOverwrite( - Scale(1 / (double) n, MulInit(&X, &hyp_tan, &temp1)), - DotMultiplyOverwrite(MulInit(&ones, Scale(A1 / (double) n, MapOverwrite(&MinusArg, A1, Sum(MapOverwrite(&Square, 0, &hyp_tan), &sum))), &temp2), &B)); + Scale(1 / (double) n, MulInit(&X, &hyp_tan, &temp1)), + DotMultiplyOverwrite(MulInit(&ones, Scale(A1 / (double) n, MapOverwrite(&MinusArg, A1, Sum(MapOverwrite(&Square, 0, &hyp_tan), &sum))), &temp2), &B)); B.PrintDebug("B"); @@ -581,5 +732,5 @@ void MakeSubMatrixByColumns(Vector column_indices, Matrix A, Matrix *A_sub) { return 0; -} + } */