|
Packit |
67cb25 |
/* tsqr.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2015 Patrick Alken
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* This module implements the sequential TSQR algorithm
|
|
Packit |
67cb25 |
* described in
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* [1] Demmel, J., Grigori, L., Hoemmen, M. F., and Langou, J.
|
|
Packit |
67cb25 |
* "Communication-optimal parallel and sequential QR and LU factorizations",
|
|
Packit |
67cb25 |
* UCB Technical Report No. UCB/EECS-2008-89, 2008.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* The algorithm operates on a tall least squares system:
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* [ A_1 ] x = [ b_1 ]
|
|
Packit |
67cb25 |
* [ A_2 ] [ b_2 ]
|
|
Packit |
67cb25 |
* [ ... ] [ ... ]
|
|
Packit |
67cb25 |
* [ A_k ] [ b_k ]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* as follows:
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* 1. Initialize
|
|
Packit |
67cb25 |
* a. [Q_1,R_1] = qr(A_1)
|
|
Packit |
67cb25 |
* b. z_1 = Q_1^T b_1
|
|
Packit |
67cb25 |
* 2. Loop i = 2:k
|
|
Packit |
67cb25 |
* a. [Q_i,R_i] = qr( [ R_{i-1} ; A_i ] )
|
|
Packit |
67cb25 |
* b. z_i = Q_i^T [ z_{i-1} ; b_i ]
|
|
Packit |
67cb25 |
* 3. Output:
|
|
Packit |
67cb25 |
* a. R = R_k
|
|
Packit |
67cb25 |
* b. Q^T b = z_k
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Step 2(a) is optimized to take advantage
|
|
Packit |
67cb25 |
* of the sparse structure of the matrix
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_vector.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_matrix.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_linalg.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_blas.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_multilarge.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_multifit.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
typedef struct
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t p; /* number of columns of LS matrix */
|
|
Packit |
67cb25 |
int init; /* QR system has been initialized */
|
|
Packit |
67cb25 |
int svd; /* SVD of R has been computed */
|
|
Packit |
67cb25 |
double normb; /* || b || for computing residual norm */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector *tau; /* Householder scalars, p-by-1 */
|
|
Packit |
67cb25 |
gsl_matrix *R; /* [ R ; A_i ], size p-by-p */
|
|
Packit |
67cb25 |
gsl_vector *QTb; /* [ Q^T b ; b_i ], size p-by-1 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_multifit_linear_workspace *multifit_workspace_p;
|
|
Packit |
67cb25 |
} tsqr_state_t;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void *tsqr_alloc(const size_t p);
|
|
Packit |
67cb25 |
static void tsqr_free(void *vstate);
|
|
Packit |
67cb25 |
static int tsqr_reset(void *vstate);
|
|
Packit |
67cb25 |
static int tsqr_accumulate(gsl_matrix * A, gsl_vector * b,
|
|
Packit |
67cb25 |
void * vstate);
|
|
Packit |
67cb25 |
static int tsqr_solve(const double lambda, gsl_vector * x,
|
|
Packit |
67cb25 |
double * rnorm, double * snorm,
|
|
Packit |
67cb25 |
void * vstate);
|
|
Packit |
67cb25 |
static int tsqr_rcond(double * rcond, void * vstate);
|
|
Packit |
67cb25 |
static int tsqr_lcurve(gsl_vector * reg_param, gsl_vector * rho,
|
|
Packit |
67cb25 |
gsl_vector * eta, void * vstate);
|
|
Packit |
67cb25 |
static int tsqr_svd(tsqr_state_t * state);
|
|
Packit |
67cb25 |
static double tsqr_householder_transform (double *v0, gsl_vector * v);
|
|
Packit |
67cb25 |
static int tsqr_householder_hv (const double tau, const gsl_vector * v, double *w0,
|
|
Packit |
67cb25 |
gsl_vector * w);
|
|
Packit |
67cb25 |
static int tsqr_householder_hm (const double tau, const gsl_vector * v, gsl_matrix * R,
|
|
Packit |
67cb25 |
gsl_matrix * A);
|
|
Packit |
67cb25 |
static int tsqr_QR_decomp (gsl_matrix * R, gsl_matrix * A, gsl_vector * tau);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
tsqr_alloc()
|
|
Packit |
67cb25 |
Allocate workspace for solving large linear least squares
|
|
Packit |
67cb25 |
problems using the TSQR approach
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: p - number of columns of LS matrix
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: pointer to workspace
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void *
|
|
Packit |
67cb25 |
tsqr_alloc(const size_t p)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
tsqr_state_t *state;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (p == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("p must be a positive integer",
|
|
Packit |
67cb25 |
GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state = calloc(1, sizeof(tsqr_state_t));
|
|
Packit |
67cb25 |
if (!state)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("failed to allocate tsqr state", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->p = p;
|
|
Packit |
67cb25 |
state->init = 0;
|
|
Packit |
67cb25 |
state->svd = 0;
|
|
Packit |
67cb25 |
state->normb = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->R = gsl_matrix_alloc(p, p);
|
|
Packit |
67cb25 |
if (state->R == NULL)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
tsqr_free(state);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("failed to allocate R matrix", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->QTb = gsl_vector_alloc(p);
|
|
Packit |
67cb25 |
if (state->QTb == NULL)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
tsqr_free(state);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("failed to allocate QTb vector", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->tau = gsl_vector_alloc(p);
|
|
Packit |
67cb25 |
if (state->tau == NULL)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
tsqr_free(state);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("failed to allocate tau vector", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->multifit_workspace_p = gsl_multifit_linear_alloc(p, p);
|
|
Packit |
67cb25 |
if (state->multifit_workspace_p == NULL)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
tsqr_free(state);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("failed to allocate multifit workspace", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return state;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void
|
|
Packit |
67cb25 |
tsqr_free(void *vstate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
tsqr_state_t *state = (tsqr_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->R)
|
|
Packit |
67cb25 |
gsl_matrix_free(state->R);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->QTb)
|
|
Packit |
67cb25 |
gsl_vector_free(state->QTb);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->tau)
|
|
Packit |
67cb25 |
gsl_vector_free(state->tau);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->multifit_workspace_p)
|
|
Packit |
67cb25 |
gsl_multifit_linear_free(state->multifit_workspace_p);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
free(state);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
tsqr_reset(void *vstate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
tsqr_state_t *state = (tsqr_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_set_zero(state->R);
|
|
Packit |
67cb25 |
gsl_vector_set_zero(state->QTb);
|
|
Packit |
67cb25 |
state->init = 0;
|
|
Packit |
67cb25 |
state->svd = 0;
|
|
Packit |
67cb25 |
state->normb = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
tsqr_accumulate()
|
|
Packit |
67cb25 |
Add a new block of rows to the QR system
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: A - new block of rows, n-by-p
|
|
Packit |
67cb25 |
b - new rhs vector n-by-1
|
|
Packit |
67cb25 |
vstate - workspace
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: success/error
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Notes:
|
|
Packit |
67cb25 |
1) On output, the upper triangular portion of state->R(1:p,1:p)
|
|
Packit |
67cb25 |
contains current R matrix
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
2) state->QTb(1:p) contains current Q^T b vector
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
3) A and b are destroyed
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
tsqr_accumulate(gsl_matrix * A, gsl_vector * b, void * vstate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
tsqr_state_t *state = (tsqr_state_t *) vstate;
|
|
Packit |
67cb25 |
const size_t n = A->size1;
|
|
Packit |
67cb25 |
const size_t p = A->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (p != state->p)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("columns of A do not match workspace", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (n != b->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("A and b have different numbers of rows", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (state->init == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
const size_t npmin = GSL_MIN(n, p);
|
|
Packit |
67cb25 |
gsl_vector_view tau = gsl_vector_subvector(state->tau, 0, npmin);
|
|
Packit |
67cb25 |
gsl_matrix_view R = gsl_matrix_submatrix(state->R, 0, 0, npmin, p);
|
|
Packit |
67cb25 |
gsl_matrix_view Av = gsl_matrix_submatrix(A, 0, 0, npmin, p);
|
|
Packit |
67cb25 |
gsl_vector_view QTb = gsl_vector_subvector(state->QTb, 0, npmin);
|
|
Packit |
67cb25 |
gsl_vector_view bv = gsl_vector_subvector(b, 0, npmin);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* this is the first matrix block A_1, compute its (dense) QR decomposition */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute QR decomposition of A */
|
|
Packit |
67cb25 |
status = gsl_linalg_QR_decomp(A, &tau.vector);
|
|
Packit |
67cb25 |
if (status)
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* store upper triangular R factor in state->R */
|
|
Packit |
67cb25 |
gsl_matrix_tricpy('U', 1, &R.matrix, &Av.matrix);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute ||b|| */
|
|
Packit |
67cb25 |
state->normb = gsl_blas_dnrm2(b);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute Q^T b and keep the first p elements */
|
|
Packit |
67cb25 |
gsl_linalg_QR_QTvec(A, &tau.vector, b);
|
|
Packit |
67cb25 |
gsl_vector_memcpy(&QTb.vector, &bv.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->init = 1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute QR decomposition of [ R_{i-1} ; A_i ], accounting for
|
|
Packit |
67cb25 |
* sparse structure */
|
|
Packit |
67cb25 |
status = tsqr_QR_decomp(state->R, A, state->tau);
|
|
Packit |
67cb25 |
if (status)
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* update ||b|| */
|
|
Packit |
67cb25 |
state->normb = gsl_hypot(state->normb, gsl_blas_dnrm2(b));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* compute Q^T [ QTb_{i - 1}; b_i ], accounting for the sparse
|
|
Packit |
67cb25 |
* structure of the Householder reflectors
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < p; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double ti = gsl_vector_get (state->tau, i);
|
|
Packit |
67cb25 |
gsl_vector_const_view h = gsl_matrix_const_column (A, i);
|
|
Packit |
67cb25 |
double *wi = gsl_vector_ptr(state->QTb, i);
|
|
Packit |
67cb25 |
tsqr_householder_hv (ti, &(h.vector), wi, b);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
tsqr_solve()
|
|
Packit |
67cb25 |
Solve the least squares system:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
chi^2 = || QTb - R x ||^2 + lambda^2 || x ||^2
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
using the SVD of R
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: lambda - regularization parameter
|
|
Packit |
67cb25 |
x - (output) solution vector p-by-1
|
|
Packit |
67cb25 |
rnorm - (output) residual norm ||b - A x||
|
|
Packit |
67cb25 |
snorm - (output) solution norm ||x||
|
|
Packit |
67cb25 |
vstate - workspace
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: success/error
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
tsqr_solve(const double lambda, gsl_vector * x,
|
|
Packit |
67cb25 |
double * rnorm, double * snorm,
|
|
Packit |
67cb25 |
void * vstate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
tsqr_state_t *state = (tsqr_state_t *) vstate;
|
|
Packit |
67cb25 |
const size_t p = x->size;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (p != state->p)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("solution vector does not match workspace", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute SVD of R if not already computed */
|
|
Packit |
67cb25 |
if (state->svd == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
status = tsqr_svd(state);
|
|
Packit |
67cb25 |
if (status)
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
status = gsl_multifit_linear_solve(lambda, state->R, state->QTb, x, rnorm, snorm,
|
|
Packit |
67cb25 |
state->multifit_workspace_p);
|
|
Packit |
67cb25 |
if (status)
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* Since we're solving a reduced square system above, we need
|
|
Packit |
67cb25 |
* to account for the full residual vector:
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* rnorm = || [ Q1^T b - R x ; Q2^T b ] ||
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* where Q1 is the thin Q factor of X, and Q2
|
|
Packit |
67cb25 |
* are the remaining columns of Q. But:
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* || Q2^T b ||^2 = ||b||^2 - ||Q1^T b||^2
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* so add this into the rnorm calculation
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double norm_Q1Tb = gsl_blas_dnrm2(state->QTb);
|
|
Packit |
67cb25 |
double ratio = norm_Q1Tb / state->normb;
|
|
Packit |
67cb25 |
double diff = 1.0 - ratio*ratio;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (diff > GSL_DBL_EPSILON)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double norm_Q2Tb = state->normb * sqrt(diff);
|
|
Packit |
67cb25 |
*rnorm = gsl_hypot(*rnorm, norm_Q2Tb);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
tsqr_lcurve()
|
|
Packit |
67cb25 |
Compute L-curve of least squares system
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: reg_param - (output) vector of regularization parameters
|
|
Packit |
67cb25 |
rho - (output) vector of residual norms
|
|
Packit |
67cb25 |
eta - (output) vector of solution norms
|
|
Packit |
67cb25 |
vstate - workspace
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: success/error
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
tsqr_lcurve(gsl_vector * reg_param, gsl_vector * rho,
|
|
Packit |
67cb25 |
gsl_vector * eta, void * vstate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
tsqr_state_t *state = (tsqr_state_t *) vstate;
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute SVD of R if not already computed */
|
|
Packit |
67cb25 |
if (state->svd == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
status = tsqr_svd(state);
|
|
Packit |
67cb25 |
if (status)
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
status = gsl_multifit_linear_lcurve(state->QTb, reg_param, rho, eta,
|
|
Packit |
67cb25 |
state->multifit_workspace_p);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* now add contribution to rnorm from Q2 factor */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double norm_Q1Tb = gsl_blas_dnrm2(state->QTb);
|
|
Packit |
67cb25 |
double ratio = norm_Q1Tb / state->normb;
|
|
Packit |
67cb25 |
double diff = 1.0 - ratio*ratio;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (diff > GSL_DBL_EPSILON)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double norm_Q2Tb = state->normb * sqrt(diff);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < rho->size; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double *rhoi = gsl_vector_ptr(rho, i);
|
|
Packit |
67cb25 |
*rhoi = gsl_hypot(*rhoi, norm_Q2Tb);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
tsqr_rcond(double * rcond, void * vstate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
tsqr_state_t *state = (tsqr_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute SVD of R if not already computed */
|
|
Packit |
67cb25 |
if (state->svd == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status = tsqr_svd(state);
|
|
Packit |
67cb25 |
if (status)
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*rcond = gsl_multifit_linear_rcond(state->multifit_workspace_p);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
tsqr_svd()
|
|
Packit |
67cb25 |
Compute the SVD of the upper triangular
|
|
Packit |
67cb25 |
R factor. This allows us to compute the upper/lower
|
|
Packit |
67cb25 |
bounds on the regularization parameter and compute
|
|
Packit |
67cb25 |
the matrix reciprocal condition number.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: state - workspace
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: success/error
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
tsqr_svd(tsqr_state_t * state)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
status = gsl_multifit_linear_svd(state->R, state->multifit_workspace_p);
|
|
Packit |
67cb25 |
if (status)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("error computing SVD of R", status);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->svd = 1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
tsqr_householder_transform()
|
|
Packit |
67cb25 |
This routine is an optimized version of
|
|
Packit |
67cb25 |
gsl_linalg_householder_transform(), designed for the QR
|
|
Packit |
67cb25 |
decomposition of M-by-N matrices of the form:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
T = [ R ]
|
|
Packit |
67cb25 |
[ A ]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where R is N-by-N upper triangular, and A is (M-N)-by-N dense.
|
|
Packit |
67cb25 |
This routine computes a householder transformation (tau,v) of a
|
|
Packit |
67cb25 |
x so that P x = [ I - tau*v*v' ] x annihilates x(1:n-1). x will
|
|
Packit |
67cb25 |
be a subcolumn of the matrix T, and so its structure will be:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
x = [ x0 ] <- 1 nonzero value for the diagonal element of R
|
|
Packit |
67cb25 |
[ 0 ] <- N - j - 1 zeros, where j is column of matrix in [0,N-1]
|
|
Packit |
67cb25 |
[ x ] <- M-N nonzero values for the dense part A
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: v0 - pointer to diagonal element of R
|
|
Packit |
67cb25 |
on input, v0 = x0;
|
|
Packit |
67cb25 |
v - on input, x vector
|
|
Packit |
67cb25 |
on output, householder vector v
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double
|
|
Packit |
67cb25 |
tsqr_householder_transform (double *v0, gsl_vector * v)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* replace v[0:M-1] with a householder vector (v[0:M-1]) and
|
|
Packit |
67cb25 |
coefficient tau that annihilate v[1:M-1] */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double alpha, beta, tau ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute xnorm = || [ 0 ; v ] ||, ignoring zero part of vector */
|
|
Packit |
67cb25 |
double xnorm = gsl_blas_dnrm2(v);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (xnorm == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return 0.0; /* tau = 0 */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
alpha = *v0;
|
|
Packit |
67cb25 |
beta = - (alpha >= 0.0 ? +1.0 : -1.0) * hypot(alpha, xnorm) ;
|
|
Packit |
67cb25 |
tau = (beta - alpha) / beta ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double s = (alpha - beta);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (fabs(s) > GSL_DBL_MIN)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_blas_dscal (1.0 / s, v);
|
|
Packit |
67cb25 |
*v0 = beta;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_blas_dscal (GSL_DBL_EPSILON / s, v);
|
|
Packit |
67cb25 |
gsl_blas_dscal (1.0 / GSL_DBL_EPSILON, v);
|
|
Packit |
67cb25 |
*v0 = beta;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return tau;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
tsqr_householder_hv()
|
|
Packit |
67cb25 |
Apply Householder reflector to a vector. The Householder
|
|
Packit |
67cb25 |
reflectors are for the QR decomposition of the matrix
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
[ R ]
|
|
Packit |
67cb25 |
[ A ]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where R is p-by-p upper triangular and A is n-by-p dense.
|
|
Packit |
67cb25 |
Therefore all relevant components of the Householder
|
|
Packit |
67cb25 |
vector are stored in the columns of A, while the components
|
|
Packit |
67cb25 |
in R are 0, except for diag(R) which are 1.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The vector w to be transformed is partitioned as
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
[ w1 ]
|
|
Packit |
67cb25 |
[ w2 ]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where w1 is p-by-1 and w2 is n-by-1. The w2 portion
|
|
Packit |
67cb25 |
of w is transformed by v, but most of w1 remains unchanged
|
|
Packit |
67cb25 |
except for the first element, w0
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: tau - Householder scalar
|
|
Packit |
67cb25 |
v - Householder vector, n-by-1
|
|
Packit |
67cb25 |
w0 - (input/output)
|
|
Packit |
67cb25 |
on input, w1(0);
|
|
Packit |
67cb25 |
on output, transformed w1(0)
|
|
Packit |
67cb25 |
w - (input/output) n-by-1
|
|
Packit |
67cb25 |
on input, vector w2;
|
|
Packit |
67cb25 |
on output, P*w2
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
tsqr_householder_hv (const double tau, const gsl_vector * v, double *w0, gsl_vector * w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* applies a householder transformation v to vector w */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (tau == 0)
|
|
Packit |
67cb25 |
return GSL_SUCCESS ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double d1, d;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute d1 = v(2:n)' w(2:n) */
|
|
Packit |
67cb25 |
gsl_blas_ddot (v, w, &d1;;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute d = v'w = w(1) + d1 since v(1) = 1 */
|
|
Packit |
67cb25 |
d = *w0 + d1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute w = w - tau (v) (v'w) */
|
|
Packit |
67cb25 |
*w0 -= tau * d;
|
|
Packit |
67cb25 |
gsl_blas_daxpy (-tau * d, v, w);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
tsqr_householder_hm()
|
|
Packit |
67cb25 |
Apply Householder reflector to a submatrix of
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
[ R ]
|
|
Packit |
67cb25 |
[ A ]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where R is p-by-p upper triangular and A is n-by-p dense.
|
|
Packit |
67cb25 |
The diagonal terms of R are already transformed by
|
|
Packit |
67cb25 |
tsqr_householder_transform(), so we just need to operate
|
|
Packit |
67cb25 |
on the submatrix A(:,i:p) as well as the superdiagonal
|
|
Packit |
67cb25 |
elements of R
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: tau - Householder scalar
|
|
Packit |
67cb25 |
v - Householder vector
|
|
Packit |
67cb25 |
R - upper triangular submatrix of R, (p-i)-by-(p-i-1)
|
|
Packit |
67cb25 |
A - dense submatrix of A, n-by-(p-i)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
tsqr_householder_hm (const double tau, const gsl_vector * v, gsl_matrix * R,
|
|
Packit |
67cb25 |
gsl_matrix * A)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* applies a householder transformation v,tau to matrix [ R ; A ] */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (tau == 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < A->size2; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double R0j = gsl_matrix_get (R, 0, j);
|
|
Packit |
67cb25 |
double wj;
|
|
Packit |
67cb25 |
gsl_vector_view A1j = gsl_matrix_column(A, j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_ddot (&A1j.vector, v, &wj;;
|
|
Packit |
67cb25 |
wj += R0j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_set (R, 0, j, R0j - tau * wj);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_daxpy (-tau * wj, v, &A1j.vector);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
tsqr_QR_decomp()
|
|
Packit |
67cb25 |
Compute the QR decomposition of the matrix
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
[ R ]
|
|
Packit |
67cb25 |
[ A ]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where R is p-by-p upper triangular and A is n-by-p dense.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: R - upper triangular p-by-p matrix
|
|
Packit |
67cb25 |
A - dense n-by-p matrix
|
|
Packit |
67cb25 |
tau - Householder scalars
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
tsqr_QR_decomp (gsl_matrix * R, gsl_matrix * A, gsl_vector * tau)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t n = A->size1;
|
|
Packit |
67cb25 |
const size_t p = R->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (R->size2 != A->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("R and A have different number of columns", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (tau->size != p)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of tau must be p", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < p; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Compute the Householder transformation to reduce the j-th
|
|
Packit |
67cb25 |
column of the matrix [ R ; A ] to a multiple of the j-th unit vector,
|
|
Packit |
67cb25 |
taking into account the sparse structure of R */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_view c = gsl_matrix_column(A, i);
|
|
Packit |
67cb25 |
double *Rii = gsl_matrix_ptr(R, i, i);
|
|
Packit |
67cb25 |
double tau_i = tsqr_householder_transform(Rii, &c.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set (tau, i, tau_i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Apply the transformation to the remaining columns and
|
|
Packit |
67cb25 |
update the norms */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i + 1 < p)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_view Rv = gsl_matrix_submatrix(R, i, i + 1, p - i, p - (i + 1));
|
|
Packit |
67cb25 |
gsl_matrix_view Av = gsl_matrix_submatrix(A, 0, i + 1, n, p - (i + 1));
|
|
Packit |
67cb25 |
tsqr_householder_hm (tau_i, &(c.vector), &(Rv.matrix), &(Av.matrix));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static const gsl_multilarge_linear_type tsqr_type =
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
"tsqr",
|
|
Packit |
67cb25 |
tsqr_alloc,
|
|
Packit |
67cb25 |
tsqr_reset,
|
|
Packit |
67cb25 |
tsqr_accumulate,
|
|
Packit |
67cb25 |
tsqr_solve,
|
|
Packit |
67cb25 |
tsqr_rcond,
|
|
Packit |
67cb25 |
tsqr_lcurve,
|
|
Packit |
67cb25 |
tsqr_free
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const gsl_multilarge_linear_type * gsl_multilarge_linear_tsqr =
|
|
Packit |
67cb25 |
&tsqr_type;
|