gmm.c

Go to the documentation of this file.

/*
 * This work is licensed under a Creative Commons 
 * Attribution-Noncommercial-Share Alike 3.0 United States License.
 * 
 *    http://creativecommons.org/licenses/by-nc-sa/3.0/us/
 * 
 * You are free:
 * 
 *    to Share - to copy, distribute, display, and perform the work
 *    to Remix - to make derivative works
 * 
 * Under the following conditions:
 * 
 *    Attribution. You must attribute the work in the manner specified by the
 *    author or licensor (but not in any way that suggests that they endorse you
 *    or your use of the work).
 * 
 *    Noncommercial. You may not use this work for commercial purposes.
 * 
 *    Share Alike. If you alter, transform, or build upon this work, you may
 *    distribute the resulting work only under the same or similar license to
 *    this one.
 * 
 * For any reuse or distribution, you must make clear to others the license
 * terms of this work. The best way to do this is by including this header.
 * 
 * Any of the above conditions can be waived if you get permission from the
 * copyright holder.
 * 
 * Apart from the remix rights granted under this license, nothing in this
 * license impairs or restricts the author's moral rights.
 */


#include <jwsc/config.h>

#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include <inttypes.h>

#include "jwsc/base/limits.h"
#include "jwsc/vector/vector.h"
#include "jwsc/vector/vector_math.h"
#include "jwsc/matrix/matrix.h"
#include "jwsc/matrix/matrix_math.h"
#include "jwsc/matblock/matblock.h"
#include "jwsc/prob/mv_pdf.h"
#include "jwsc/stat/stat.h"
#include "jwsc/stat/kmeans.h"
#include "jwsc/stat/gmm.h"


static void ensure_non_singular_covariance_d(Matrix_d** m, double min)
{
    uint32_t n, N;

    Matrix_d* v  = NULL;
    Vector_d* d  = NULL;
    Matrix_d* dd = NULL;
    Error*    e;

    N = (*m)->num_rows;

    /*
    for (n = 0; n < N; n++)
    {
        if ((*m)->elts[ n ][ n ] < min)
        {
            (*m)->elts[ n ][ n ] = min;
        }
    }
    */

    assert((e = solve_real_symmetric_matrix_eigenproblem_d(&v, &d, *m)) == NULL);
    create_zero_matrix_d(&dd, N, N);

    for (n = 0; n < N; n++)
    {
        if (d->elts[ n ] < min)
        {
            dd->elts[ n ][ n ] = min;
        }
        else
        {
            dd->elts[ n ][ n ] = d->elts[ n ];
        }
    }

    multiply_matrices_d(m, v, dd);
    transpose_matrix_d(&v, v);
    multiply_matrices_d(m, *m, v);

    free_matrix_d(v);
    free_vector_d(d);
    free_matrix_d(dd);
}


static void gmm_e_step_d
(
    const Matrix_d*   mu,
    const Matblock_d* sigma,
    const Vector_d*   pi,
    Matrix_d*         gamma,
    const Matrix_d*   x
)
{
    uint32_t n, N;
    uint32_t k, K;
    double   C_n;

    Vector_d* p_k     = 0;
    Vector_d* mu_k    = 0;
    Matrix_d* sigma_k = 0;
    Error*    e;

    N = x->num_rows;
    K = pi->num_elts;

    create_zero_matrix_d(&gamma, N, K);

    for (k = 0; k < K; k++)
    {
        copy_vector_from_matrix_d(&mu_k, mu, MATRIX_ROW_VECTOR, k);
        copy_matrix_from_matblock_d(&sigma_k, sigma, MATBLOCK_MAT_MATRIX, k);

        assert((e = set_mv_gaussian_pdf_d(&p_k, mu_k, sigma_k, x)) == 0);

        copy_vector_into_matrix_d(&gamma, gamma, p_k, MATRIX_COL_VECTOR, k);
    }

    for (n = 0; n < N; n++)
    {
        C_n = 0;

        for (k = 0; k < K; k++)
        {
            gamma->elts[ n ][ k ] *= pi->elts[ k ];

            C_n += gamma->elts[ n ][ k ];
        }

        assert(C_n > 0);

        for (k = 0; k < K; k++)
        {
            gamma->elts[ n ][ k ] /= C_n;
        }
    }

    free_vector_d(p_k);
    free_vector_d(mu_k);
    free_matrix_d(sigma_k);
}


static void gmm_m_step_d
(
    Matrix_d*       mu,
    Matblock_d*     sigma,
    Vector_d*       pi,
    const Matrix_d* gamma,
    const Matrix_d* x
)
{
    uint32_t k, K;
    uint32_t n, N;
    uint32_t d, dd, D;
    double   s_1, s_2;

    Vector_d* mu_k    = NULL;
    Matrix_d* sigma_k = NULL;

    K = pi->num_elts;
    N = x->num_rows;
    D = x->num_cols;

    create_zero_vector_d(&pi, K);

    // mu
    create_matrix_d(&mu, K, D);
    for (k = 0; k < K; k++)
    {
        create_zero_vector_d(&mu_k, D);

        for (n = 0; n < N; n++)
        {
            for (d = 0; d < D; d++)
            {
                mu_k->elts[ d ] += gamma->elts[ n ][ k ]*x->elts[ n ][ d ];
            }

            pi->elts[ k ] += gamma->elts[ n ][ k ];
        }

        assert(pi->elts[ k ] > 0);

        multiply_vector_by_scalar_d(&mu_k, mu_k, 1.0 / pi->elts[ k ]);
        copy_vector_into_matrix_d(&mu, mu, mu_k, MATRIX_ROW_VECTOR, k);
    }

    // sigma
    create_matblock_d(&sigma, K, D, D);
    for (k = 0; k < K; k++)
    {
        create_zero_matrix_d(&sigma_k, D, D);

        for (n = 0; n < N; n++)
        {
            for (d = 0; d < D; d++)
            {
                s_1 = x->elts[ n ][ d ] - mu->elts[ k ][ d ];

                for (dd = 0; dd < D; dd++)
                {
                    s_2 = x->elts[ n ][ dd ] - mu->elts[ k ][ dd ];

                    sigma_k->elts[ d ][ dd ] += gamma->elts[ n ][ k ] * s_1*s_2;
                }
            }
        }

        multiply_matrix_by_scalar_d(&sigma_k, sigma_k, 1.0 / pi->elts[ k ]);
        ensure_non_singular_covariance_d(&sigma_k, 1.0e-1);
        copy_matrix_into_matblock_d(&sigma, sigma, sigma_k, 
                MATBLOCK_MAT_MATRIX, k);
    }

    // pi
    multiply_vector_by_scalar_d(&pi, pi, 1.0 / N);

    free_vector_d(mu_k);
    free_matrix_d(sigma_k);
}


double gmm_log_likelihood_d
(
    const Matrix_d*   mu,
    const Matblock_d* sigma,
    const Vector_d*   pi,
    const Matrix_d*   x
)
{
    uint32_t n, N;
    uint32_t k, K;
    double   ll;

    Matrix_d* p       = 0;
    Vector_d* p_k     = 0;
    Vector_d* q       = 0;
    Vector_d* mu_k    = 0;
    Matrix_d* sigma_k = 0;
    Error*    e;

    N = x->num_rows;
    K = pi->num_elts;

    create_matrix_d(&p, N, K);

    for (k = 0; k < K; k++)
    {
        copy_vector_from_matrix_d(&mu_k, mu, MATRIX_ROW_VECTOR, k);
        copy_matrix_from_matblock_d(&sigma_k, sigma, MATBLOCK_MAT_MATRIX, k);

        assert((e = set_mv_gaussian_pdf_d(&p_k, mu_k, sigma_k, x)) == 0);

        copy_vector_into_matrix_d(&p, p, p_k, MATRIX_COL_VECTOR, k);
    }

    multiply_matrix_and_vector_d(&q, p, pi);
    ll = 0;

    for (n = 0; n < N; n++)
    {
        ll += log((q->elts[ n ] < JWSC_MIN_LOG_ARG) 
                ? JWSC_MIN_LOG_ARG : q->elts[ n ]);
    }

#if defined(JWSC_HAVE_ISNAN) && defined(JWSC_HAVE_ISINF)
    assert(!isnan(ll) && !isinf(ll));
#endif

    free_matrix_d(p);
    free_vector_d(p_k);
    free_vector_d(q);
    free_vector_d(mu_k);
    free_matrix_d(sigma_k);

    return ll;
}


void gmm_kmeans_init_d
(
    Matrix_d**      mu_out,
    Matblock_d**    sigma_out,
    Vector_d**      pi_out,
    const Matrix_d* x,
    uint32_t         K
)
{
    uint32_t d, D;
    uint32_t n, nn, N;
    uint32_t k;

    Vector_i32* members = NULL;
    Vector_u32* counts  = NULL;
    Matrix_d*   x_k     = NULL;
    Matrix_d*   sigma_k = NULL;

    D = x->num_cols;
    N = x->num_rows;

    init_kmeans_randomly_d(mu_out, x, K);
    kmeans_d(mu_out, &members, x, K, 2);

    create_matblock_d(sigma_out, K, D, D);
    create_zero_vector_u32(&counts, K);

    for (n = 0; n < N; n++)
    {
        k = members->elts[ n ];
        counts->elts[ k ]++;
    }
    for (k = 0; k < K; k++)
    {
        create_matrix_d(&x_k, counts->elts[ k ], D);

        nn = 0;
        for (n = 0; n < N; n++)
        {
            if (k == members->elts[ n ])
            {
                for (d = 0; d < D; d++)
                {
                    x_k->elts[ nn ][ d ] = x->elts[ n ][ d ];
                }
                nn++;
            }
        }

        mv_sample_covariance_d(&sigma_k, x_k);
        ensure_non_singular_covariance_d(&sigma_k, 1.0e-1);
        copy_matrix_into_matblock_d(sigma_out, *sigma_out, sigma_k,
                MATBLOCK_MAT_MATRIX, k);
    }

    create_vector_d(pi_out, K);
    for (k = 0; k < K; k++)
    {
        (*pi_out)->elts[ k ] = counts->elts[ k ] / (double)N;
    }

    free_vector_i32(members);
    free_vector_u32(counts);
    free_matrix_d(x_k);
    free_matrix_d(sigma_k);
}


void gmm_stats_init_d
(
    Matrix_d**      mu_out,
    Matblock_d**    sigma_out,
    Vector_d**      pi_out,
    const Matrix_d* x,
    uint32_t        K
)
{
    uint32_t k;
    uint32_t D;

    Vector_d* mu_0    = 0;
    Vector_d* mu_k    = 0;
    Matrix_d* sigma_k = 0;

    D = x->num_cols;

    if (*pi_out == NULL)
    {
        create_init_vector_d(pi_out, K, 1.0 / (double)K);
    }
    else
    {
        assert((*pi_out)->num_elts == K);
    }

    if (*sigma_out == NULL)
    {
        create_matblock_d(sigma_out, K, D, D);
        mv_sample_covariance_d(&sigma_k, x);
        ensure_non_singular_covariance_d(&sigma_k, 1.0e-1);

        for (k = 0; k < K; k++)
        {
            copy_matrix_into_matblock_d(sigma_out, *sigma_out, sigma_k,
                    MATBLOCK_MAT_MATRIX, k);
        }
    }
    else
    {
        assert((*sigma_out)->num_mats == K && (*sigma_out)->num_rows == D 
                                           && (*sigma_out)->num_cols == D);

        for (k = 0; k < K; k++)
        {
            copy_matrix_from_matblock_d(&sigma_k, *sigma_out, 
                    MATBLOCK_MAT_MATRIX, k);
            ensure_non_singular_covariance_d(&sigma_k, 1.0e-1);
            copy_matrix_into_matblock_d(sigma_out, *sigma_out, sigma_k,
                    MATBLOCK_MAT_MATRIX, k);
        }
    }

    if (*mu_out == NULL)
    {
        create_matrix_d(mu_out, K, D);
        mv_sample_mean_d(&mu_0, x);
        copy_vector_into_matrix_d(mu_out, *mu_out, mu_0, MATRIX_ROW_VECTOR, 0);

        for (k = 1; k < K; k++)
        {
            sample_mv_gaussian_pdf_d(&mu_k, mu_0, sigma_k);
            copy_vector_into_matrix_d(mu_out, *mu_out, mu_k, 
                    MATRIX_ROW_VECTOR, k);
        }
    }
    else
    {
        assert((*mu_out)->num_rows == K && (*mu_out)->num_cols == D);
    }

    free_vector_d(mu_0);
    free_vector_d(mu_k);
    free_matrix_d(sigma_k);
}


void train_gmm_d
(
    Matrix_d**      mu_out,
    Matblock_d**    sigma_out,
    Vector_d**      pi_out,
    Matrix_d**      gamma_out,
    const Matrix_d* x,
    uint32_t        K
)
{
    double   ll, ll_prev;

    Matrix_d*   mu;
    Matblock_d* sigma;
    Vector_d*   pi;
    Matrix_d*   gamma = 0;

    gmm_kmeans_init_d(mu_out, sigma_out, pi_out, x, K);
    mu    = *mu_out;
    sigma = *sigma_out;
    pi    = *pi_out;

    if (gamma_out)
    {
        create_matrix_d(gamma_out, x->num_rows, pi->num_elts);
        gamma = *gamma_out;
    }
    else
    {
        create_matrix_d(&gamma, x->num_rows, pi->num_elts);
    }

    ll_prev = 0;
    ll      = log(0);

    while (fabs(ll - ll_prev) > 1.0e-10)
    {
        gmm_e_step_d(mu, sigma, pi, gamma, x);
        gmm_m_step_d(mu, sigma, pi, gamma, x);

        ll_prev = ll;
        ll = gmm_log_likelihood_d(mu, sigma, pi, x);
    }

    if (!gamma_out)
    {
        free_matrix_d(gamma);
    }
}