SuiteSparseQR_Cへの道 その９ 平滑化拘束条件付きBスプライン + ABIC

前回までで連立方程式を組み立てる準備が整いました。
今回はそれを解くルーチンを用意しましょう。

今回することは、観測値ベクトル y、行列 z 及び解ベクトル x について

連立方程式  y = z * x

をQR分解を使って解きます。
但し z は疎行列で。

また、ABICを計算する為に z のランクと行列式の log も計算します。

QR分解には関数

long
SuiteSparseQR_C_QR (int ordering, int tol, size_t econ, cholmod_sparse *z,
                             cholmod_sparse **q, cholmod_sparse **r, long **e,
                             cholmod_common *cc)

を使います。
この関数の返り値が z の rank になります。

やる事をOctaveやMatlab風に書くと

[q, r] = qr(z);  (QR分解)
qty=q' * y;   (qty = Q^T * y)
x = r \ qty; (backslash)

で解 x を求めます。

backslash の部分は関数

cholmod_dense *
SuiteSparseQR_C_backslash_default (cholmod_sparse *r, cholmod_dense *qty,
                                   cholmod_common *cc)

で求めます。
これの返り値が解ベクトルになります。
default じゃなく、order や tol を決めても良いんでしょうけどね。

また、計算の途中で行列 z のrankと log | det | も求めます。

#include <SuiteSparseQR_C.h>

enum {
    CholmodNoTrans,
    CholmodTrans,
    CholmodNTransEntries
};

enum {
    Cholmod0norm,
    Cholmod1norm,
    Cholmod2norm,
    CholmodNnormEntries
};

double    zero[] = {0., 0.};
double    one[] = {1., 0.};
double    minusone[] = {-1., 0.};

/* 行列式の log を計算 */
double
_log_det (long rank, cholmod_sparse *r)
{
    int     j, p;
    long    *ri = (long *) r->i;
    long    *rp = (long *) r->p;
    double  *rx = (double *) r->x;
    double  log_det = 0.;
    for (j = 0; j < (int) rank; j++) {
        for (p = (int) rp[j]; p < (int) rp[j + 1]; p++) {
            if ((int) ri[p] == j) {
                log_det += log (fabs (rx[p]));
                break;
            }
        }
    }
    return log_det;
}

/* 並べ替え */
void
_inv_permutate (long *e, cholmod_dense *x, cholmod_common *cc)
{
    int              i;
    cholmod_dense    *tmp = cholmod_l_copy_dense (x, cc);
    for (i = 0; i < (int) x->nzmax; i++)
        ((double *) x->x)[(size_t) e[i]] = ((double *) tmp->x)[i];
    cholmod_l_free_dense (&tmp, cc);
    return;
}

/* 方程式 y = z * x をQR分解を用いて解く */
cholmod_dense *
_solve_simeq (cholmod_sparse *z, cholmod_dense *y, long *rank, double *log_det, cholmod_common *cc)
{
    long              *e = NULL;
    cholmod_sparse    *q;
    cholmod_sparse    *r;
    cholmod_dense     *qty;
    cholmod_dense     *x;

    /* QR分解 : [q,r,p] = qr (z, 'econ'); */
    *rank = SuiteSparseQR_C_QR (SPQR_ORDERING_DEFAULT, SPQR_DEFAULT_TOL, 0, z, &q, &r, &e, cc);

    /* zの行列式 (= r の行列式) を計算 */
    *log_det = _log_det (*rank, r);

    /* qty = q' * y */
    qty = cholmod_l_allocate_dense (q->ncol, 1, q->ncol, CHOLMOD_REAL, cc);
    cholmod_l_sdmult (q, CholmodTrans, one, zero, y, qty, cc);
    cholmod_l_free_sparse (&q, cc);

    /* x = r \ (q' * y) */
    x = SuiteSparseQR_C_backslash_default (r, qty, cc);
    cholmod_l_free_sparse (&r, cc);
    cholmod_l_free_dense (&qty, cc);

    if (e) {
        _inv_permutate (e, x, cc);
        free (e);
    }

    return x;
}

/* sigma2 = |z * e - y|^2 / N */
double
_sigma2 (cholmod_dense *y, cholmod_sparse *z, cholmod_dense *x, cholmod_common *cc)
{
    double           res;
    cholmod_dense    *r = cholmod_l_copy_dense (y, cc);
    cholmod_l_sdmult (z, CholmodNoTrans, one, minusone, x, r, cc);
    res = pow (cholmod_l_norm_dense (r, Cholmod2norm, cc), 2.) / (double) r->nzmax;
    cholmod_l_free_dense (&r, cc);
    return res;
}

最初の関数

double _log_det (long rank, cholmod_sparse *r)

では、r の対角成分の絶対値の log を、r の rank 行目まで和を取っています。

void _inv_permutate (long *e, cholmod_dense *x, cholmod_common *cc)

では、QR分解で得られた解 cholmod_dense *x を、permutation である long *e に従って並べ直しています。

続く関数

cholmod_dense *
_solve_simeq (cholmod_sparse *z, cholmod_dense *y, long *rank, double *log_det,
                  cholmod_common *cc)

で、連立方程式 y = z * x を解いて、解 cholmod_dense *x を返します。
また z のランクと行列式を、３、４番目の引数 long *rank, double *log_det に格納します。

最後の

double
_sigma2 (cholmod_dense *y, cholmod_sparse *z, cholmod_dense *x, cholmod_common *cc)

は、現時点ではどこでも呼ばれませんが、計算値 (z * x) と観測値 (y) の残差の分散を求めています。これはいずれ ABIC を計算する際に必要になります。

まぁ、こんなもんでどうでしょう？
問題あるかなぁ・・・（汗