tests : add SVD experiments

diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index a20cf2859441a1205034add1bb142c4e33e2bbaf..856cc3ed8572ef6d933b5f1dee9fb8d812e1edd4 100644 (file)
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -104,3 +104,15 @@ set(TEST_TARGET test3)
 add_executable(${TEST_TARGET} ${TEST_TARGET}.c)
 target_link_libraries(${TEST_TARGET} PRIVATE ggml)
 add_test(NAME ${TEST_TARGET} COMMAND $<TARGET_FILE:${TEST_TARGET}>)
+
+#
+# test-svd0 (arm)
+
+if (${CMAKE_SYSTEM_PROCESSOR} MATCHES "arm" AND NOT GGML_NO_ACCELERATE)
+    set(TEST_TARGET test-svd0)
+    add_executable(${TEST_TARGET} ${TEST_TARGET}.c)
+    target_link_libraries(${TEST_TARGET} PRIVATE ggml ${GGML_EXTRA_LIBS})
+    target_compile_options(${TEST_TARGET} PRIVATE ${GGML_EXTRA_FLAGS})
+    add_test(NAME ${TEST_TARGET} COMMAND $<TARGET_FILE:${TEST_TARGET}>)
+endif()
+

diff --git a/tests/test-svd0.c b/tests/test-svd0.c
new file mode 100644 (file)
index 0000000..10539ba
--- /dev/null
+++ b/tests/test-svd0.c
@@ -0,0 +1,218 @@
+// SVD dimensionality reduction
+
+#include <float.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <assert.h>
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
+#include <math.h>
+
+#include <sys/time.h>
+
+#ifdef GGML_USE_ACCELERATE
+#include <Accelerate/Accelerate.h>
+#endif
+
+float frand() {
+    return (float) rand() / (float) RAND_MAX;
+}
+
+//int sgesvd_(char *__jobu, char *__jobvt, __CLPK_integer *__m,
+//        __CLPK_integer *__n, __CLPK_real *__a, __CLPK_integer *__lda,
+//        __CLPK_real *__s, __CLPK_real *__u, __CLPK_integer *__ldu,
+//        __CLPK_real *__vt, __CLPK_integer *__ldvt, __CLPK_real *__work,
+//        __CLPK_integer *__lwork,
+//        __CLPK_integer *__info)
+
+int main(int argc, const char ** argv) {
+    int m = 10;
+    int n = 5;
+
+    float * A  = (float *) malloc(n * m * sizeof(float));
+    float * A0 = (float *) malloc(n * m * sizeof(float));
+
+    for (int i = 0; i < n; ++i) {
+        for (int j = 0; j < m; ++j) {
+            A[i * m + j] = (float) (10.0f*(i + 1) + 1.0f * frand());
+            //A[i * m + j] = (float) (10.0f*(i%2 + 1) + 0.1f * frand());
+            //if (i == 2) {
+            //    A[i * m + j] += 20*frand();
+            //}
+            if ((i == 1 || i == 3) && j > m/2) {
+                A[i * m + j] = -A[i * m + j];
+            }
+        }
+    }
+
+    // average vector
+    //float * M = (float *) malloc(m * sizeof(float));
+
+    //{
+    //    for (int j = 0; j < m; ++j) {
+    //        M[j] = 0.0f;
+    //    }
+    //    for (int i = 0; i < n; ++i) {
+    //        for (int j = 0; j < m; ++j) {
+    //            M[j] += A[i * m + j];
+    //        }
+    //    }
+    //    for (int j = 0; j < m; ++j) {
+    //        M[j] /= (float) n;
+    //    }
+    //}
+
+    //// subtract average vector
+    //for (int i = 0; i < n; ++i) {
+    //    for (int j = 0; j < m; ++j) {
+    //        A[i * m + j] -= M[j];
+    //    }
+    //}
+
+    memcpy(A0, A, n * m * sizeof(float));
+
+    // print A
+    printf("A:\n");
+    for (int i = 0; i < n; ++i) {
+        printf("col %d : ", i);
+        for (int j = 0; j < m; ++j) {
+            printf("%9.5f ", A[i * m + j]);
+        }
+        printf("\n");
+    }
+    printf("\n");
+
+    // SVD
+    // A = U * S * V^T
+
+    float * U = (float *) malloc(n * m * sizeof(float));
+    float * S = (float *) malloc(n * sizeof(float));
+    float * V = (float *) malloc(n * n * sizeof(float));
+
+    int lda = m;
+    int ldu = m;
+    int ldvt = n;
+
+    float work_size;
+    int lwork = -1;
+    int info = 0;
+
+    sgesvd_("S", "S", &m, &n, A, &lda, S, U, &ldu, V, &ldvt, &work_size, &lwork, &info);
+
+    lwork = (int) work_size;
+
+    printf("work_size = %f, info = %d, lwork = %d\n", work_size, info, lwork);
+
+    float * work = (float *) malloc(lwork * sizeof(float));
+
+    sgesvd_("S", "S", &m, &n, A, &lda, S, U, &ldu, V, &ldvt, work, &lwork, &info);
+
+    // print U
+    printf("U:\n");
+    for (int i = 0; i < n; ++i) {
+        printf("col %d : ", i);
+        for (int j = 0; j < m; ++j) {
+            printf("%9.5f ", U[i * m + j]);
+        }
+        printf("\n");
+    }
+    printf("\n");
+
+    // normalize S
+    {
+        double sum = 0.0;
+        for (int i = 0; i < n; ++i) {
+            sum += S[i];
+        }
+        sum *= sqrt((double) m);
+        for (int i = 0; i < n; ++i) {
+            S[i] /= sum;
+        }
+    }
+
+    // print S
+    printf("S:\n");
+    for (int i = 0; i < n; ++i) {
+        printf("- %d = %9.5f\n", i, S[i]);
+    }
+    printf("\n");
+
+    // print V
+    printf("V:\n");
+    for (int i = 0; i < n; ++i) {
+        printf("col %d : ", i);
+        for (int j = 0; j < n; ++j) {
+            printf("%9.5f ", V[i * n + j]);
+        }
+        printf("\n");
+    }
+    printf("\n");
+
+    // print A
+    printf("A:\n");
+    for (int i = 0; i < n; ++i) {
+        printf("col %d : ", i);
+        for (int j = 0; j < m; ++j) {
+            printf("%9.5f ", A[i * m + j]);
+        }
+        printf("\n");
+    }
+    printf("\n");
+
+    // compute singular vectors in U
+    for (int i = 0; i < n; ++i) {
+        for (int j = 0; j < m; ++j) {
+            U[i * m + j] *= S[i];
+        }
+    }
+
+    // normalize U
+    for (int i = 0; i < n; ++i) {
+        double sum = 0.0;
+        for (int j = 0; j < m; ++j) {
+            sum += U[i * m + j] * U[i * m + j];
+        }
+        sum = sqrt(sum);
+        for (int j = 0; j < m; ++j) {
+            U[i * m + j] /= sum*sqrt((double) m);
+        }
+    }
+
+    // print U
+    printf("U:\n");
+    for (int i = 0; i < n; ++i) {
+        printf("col %d : ", i);
+        for (int j = 0; j < m; ++j) {
+            printf("%9.5f ", U[i * m + j]);
+        }
+        printf("\n");
+    }
+    printf("\n");
+
+
+    // project A0 onto U
+    float * A1 = (float *) malloc(n * n * sizeof(float));
+
+    for (int i = 0; i < n; ++i) {
+        for (int j = 0; j < n; ++j) {
+            A1[i * n + j] = 0.0f;
+            for (int k = 0; k < m; ++k) {
+                A1[i * n + j] += A0[i * m + k] * U[j * m + k];
+            }
+        }
+    }
+
+    // print A1
+    printf("A1:\n");
+    for (int i = 0; i < n; ++i) {
+        printf("col %d : ", i);
+        for (int j = 0; j < n; ++j) {
+            printf("%9.5f ", A1[i * n + j]);
+        }
+        printf("\n");
+    }
+    printf("\n");
+
+    return 0;
+}