gctl_optimization/lib/optimization/svd.cpp

/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#include "svd.h"

double my_norm(gctl::array<double> &arr, double eps)
{
    double norm_val = 0.0;
    for (size_t i = 0; i < arr.size(); i++)
    {
        norm_val += arr[i]*arr[i];
    }
    norm_val = sqrt(norm_val);
    
    if (norm_val < eps) return 0.0;

    for (size_t i = 0; i < arr.size(); i++)
    {
        arr[i] /= norm_val;
    }
    return norm_val;
}
    
gctl::svd::svd()
{
    reset();
}

gctl::svd::svd(const matrix<double> &src_mat) : svd()
{
    decompose(src_mat);
}

void gctl::svd::reset()
{
    maxi_iteration = 1000;
    K = 0;
    epsilon = 1e-8;
    U.clear();
    V.clear();
    S.clear();
    return;
}

void gctl::svd::set_singular_number(int k)
{
    if (k <= 0)
    {
        throw invalid_argument("Invalid singular number. From gctl::svd::set_singular_number(...)");
    }

    K = k;
    return;
}

void gctl::svd::set_iteration(int t)
{
    if (t <= 0)
    {
        throw invalid_argument("Invalid singular number. From gctl::svd::set_iteration(...)");
    }

    maxi_iteration = t;
    return;
}

void gctl::svd::set_epsilon(double e)
{
    if (e <= 0)
    {
        throw invalid_argument("Invalid singular number. From gctl::svd::set_epsilon(...)");
    }

    epsilon = e;
    return;
}

void gctl::svd::decompose(const matrix<double> &src_mat)
{
    int M = src_mat.row_size();
    int N = src_mat.col_size();
    if (K == 0) K = N;

    S.resize(K, 0.0);
    U.resize(K, M, 0.0);
    V.resize(K, N, 0.0);

    srand(time(0));
    array<double> left_vector(M), next_left_vector(M);
    array<double> right_vector(N), next_right_vector(N);
    array<double> U_tmp(M), V_tmp(N);

    double diff, r, d;
    for(int col=0;col<K;col++)
    {
        diff = 1;
        r = -1;
        while(1)
        {
            for(int i=0;i<M;i++)
                left_vector[i]= (double) rand() / RAND_MAX;
            if(my_norm(left_vector, epsilon) > epsilon)
                break;
        }

        for(int iter=0; diff >= epsilon && iter < maxi_iteration; iter++)
        {
            next_left_vector.assign_all(0.0);
            next_right_vector.assign_all(0.0);

            for(int i=0;i<M;i++)
                for(int j=0;j<N;j++)
                    next_right_vector[j]+=left_vector[i]*src_mat[i][j];

            r = my_norm(next_right_vector, epsilon);
            if(r<epsilon) break;

            for(int i=0;i<col;i++)
            {
                for (int j = 0; j < N; j++)
                {
                    V_tmp[j] = V[i][j];
                }
                V_tmp.orth(next_right_vector);
            }
            my_norm(next_right_vector, epsilon);

            for(int i=0;i<M;i++)
                for(int j=0;j<N;j++)
                    next_left_vector[i]+=next_right_vector[j]*src_mat[i][j];
            r = my_norm(next_left_vector, epsilon);
            if(r<epsilon) break;

            for(int i=0;i<col;i++)
            {
                for (int j = 0; j < M; j++)
                {
                    U_tmp[j] = U[i][j];
                }
                U_tmp.orth(next_left_vector);
            }
            my_norm(next_left_vector, epsilon);

            diff=0;
            for(int i=0;i<M;i++)
            {
                d=next_left_vector[i]-left_vector[i];
                diff+=d*d;
            }

            for (int i = 0; i < M; i++)
            {
                left_vector[i] = next_left_vector[i];
            }

            for (int i = 0; i < N; i++)
            {
                right_vector[i] = next_right_vector[i];
            }
        }

        if(r>=epsilon)
        {
            S[col]=r;
            for (int i = 0; i < M; i++)
            {
                U[col][i] = left_vector[i];
            }

            for (int i = 0; i < N; i++)
            {
                V[col][i] = right_vector[i];
            }
        }
        else break;
    }
    return;
}