arpack_svd

more...

Computes some of the singular values and singular vectors of a real rectangular matrix A = UΣVT.

Note: Function arpack_svd is available in double precision only.

Synopsis

#include <imsl.h>

double *imsl_d_arpack_svd (void fcn(), int m, int n, int nsv, …, 0)

Required Arguments

void fcn (int m, int n, double x[], int task, double y[]) (Input)

User-supplied function to return matrix-vector operations.

int m (Input)

The number of rows of matrix A.

int n (Input)

The number of columns of matrix A.

double x[] (Input)

An array of size n or m containing the vector x of the matrix-vector product Ax or AT x.

int task (Input)

An enumeration type which specifies the operation to be performed. Variable task is an enumerated integer value associated with enum type Imsl_arpack_task.

The following values are possible:

Tasks IMSL_ARPACK_A_X and IMSL_ARPACK_TRANS_A_X always have to be defined within the function.

double y[] (Output)

An array of size m or n containing the result of the matrix-vector operation Ax or AT x.

int m (Input)

The number of rows of matrix A.

int n (Input)

The number of columns of matrix A.

int nsv (Input)

The number of singular values to be computed. It is required that 1 <= nsv <= min(m,n).

Return Value

A pointer to the nsv singular values of the matrix A = UΣVT in decreasing order. To release this space, use imsl_free. If no value can be computed, then NULL is returned.

Synopsis with Optional Arguments

#include <imsl.h>

double *imsl_d_arpack_svd (void fcn(), int m, int n, int nsv,

IMSL_XGUESS, double xguess[],

IMSL_ITMAX, int itmax,

IMSL_TOLERANCE, double tol,

IMSL_SINGVAL_LOCATION, Imsl_arpack_eigval_location singval_loc,

IMSL_MATRIX_TYPE, Imsl_arpack_matrix_type matrix_type,

IMSL_NUM_LANCZOS_VECTORS, int ncv,

IMSL_NUM_ACCURATE_SINGVALS, int *n_acc,

IMSL_U, double **p_u,

IMSL_U_USER, double u[],

IMSL_U_COL_DIM, int u_col_dim,

IMSL_V, double **p_v,

IMSL_V_USER, double v[],

IMSL_V_COL_DIM, int v_col_dim,

IMSL_RETURN_USER, double s[],

IMSL_FCN_W_DATA, void fcn(), void *data,

0)

Optional Arguments

IMSL_XGUESS, double xguess[] (Input)
A non-zero vector containing the starting vector for the implicitly restarted Arnoldi method. The size of xguess is min(m,n) if matrix_type = IMSL_ARPACK_NORMAL and m + n if matrix_type = IMSL_ARPACK_EXPANDED.

By default, a random starting vector is computed internally.

IMSL_ITMAX, int itmax (Input)
The maximum number of Arnoldi iterations.
Default: itmax = 1000.

IMSL_TOLERANCE, double tol (Input)
Tolerance value used in the criterion for the acceptability of the relative accuracy of the Ritz values. A value tol<= 0 is reset to the default value.
Default: tol = imsl_d_machine(3).

IMSL_SINGVAL_LOCATION, Imsl_arpack_eigval_location singval_loc (Input)
An enumeration type which specifies the location of the singular values to compute.

Default: singval_loc = IMSL_ARPACK_LARGEST_ALGEBRAIC.

IMSL_MATRIX_TYPE, Imsl_arpack_matrix_type matrix_type (Input)
An enumeration type indicating which matrix type is used for the singular value computation.

Default: matrix_type = IMSL_ARPACK_NORMAL.

IMSL_NUM_LANCZOS_VECTORS, int ncv (Input)
The number of Lanczos vectors generated in each iteration of the Arnoldi method. It is required that nsv + 1 <= ncv <= up, where up = min(m,n) if matrix_type = IMSL_ARPACK_NORMAL and up = m + n if matrix_type = IMSL_ARPACK_EXPANDED. A value ncv >= min(2*nsv, up) is recommended.

Default: ncv = min(2*nsv, up).

IMSL_NUM_ACCURATE_SINGVALS, int *n_acc (Output)
The number of singular values that the algorithm was able to compute accurately. This number can be smaller than nsv.

IMSL_U, double **p_u (Output)
The address of a pointer to an array of size n × nsv containing the left singular vectors of A in the first n_acc columns. Typically, double * p_u is declared, and &p_u is used as an argument.

In general, the computability of a singular vector depends on whether the associated singular value is larger than zero or zero. The following table shows which singular vectors can be computed:

IMSL_U_USER, double u[](Output)
A user-defined array of size m×nsv containing the left singular vectors of A in the first n_acc columns.

IMSL_U_COL_DIM, int u_col_dim (Input)

The column dimension of the array u containing the left singular vectors.
Default: u_col_dim = nsv.

IMSL_V, double **p_v (Output)
The address of a pointer to an array of size n× nsv containing the right singular vectors of A in the first n_acc columns. On return, the necessary space is allocated by imsl_d_arpack_svd. Typically, double *p_v is declared, and &p_v is used as an argument.

In general, the computability of a singular vector depends on whether the associated singular value is larger than zero or zero.

The following table shows which singular vectors can be computed:

IMSL_V_USER, double v[] (Output)
A user-defined array of size n×nsv containing the right singular vectors of A in the first n_acc columns.

IMSL_V_COL_DIM, int v_col_dim (Input)

The column dimension of the array v containing the right singular vectors.
Default: v_col_dim = nsv.

IMSL_RETURN_USER, double s[] (Output)
An array of size nsv containing the accurately computed singular values in decreasing order in the first n_acc locations.

IMSL_FCN_W_DATA, void fcn (int m, int n, double x[], int task, double y[], void *data), void *data (Input/Output)
User-supplied function to return matrix-vector operations, which also accepts a pointer to data that is supplied by the user. data is a pointer to the data to be passed to the user-supplied function. See Passing Data to User-Supplied Functions in the introduction to this manual for more details.

Note: The possibility to supply user-data via IMSL_FCN_W_DATA is an important feature of arpack_svd. It allows the user to transfer problem-specific data to the algorithm without the need to define global data. See the documentation example 2 for a use-case.

Description

Function imsl_d_arpack_svd computes a partial singular value decomposition of a general real-valued m×n matrix A. The singular value decomposition (SVD) of A is defined as

A = UΣVT,

where , and

The σi are called the singular values of A. Matrices U∈Rm×m and V∈Rn×n are orthogonal, and setting

the ui and vi are, respectively, the left and right singular vectors associated with σi , i=1,…, min (m,n).

Function imsl_d_arpack_svd offers two simple transformations that allow to convert results about eigenvalues and eigenvectors of symmetric matrices into results about singular values and singular vectors of arbitrary real matrices.

In the normal approach, selectable by setting matrix_type = IMSL_ARPACK_NORMAL, function imsl_d_arpack_svd iteratively computes selected eigenvalues and eigenvectors of the symmetric matrices

ATA, if m≥n,

AAT if m<n,

by calling function imsl_d_arpack_symmetric internally. Since

the non-negative square roots of the eigenvalues of ATA correspond to the singular values of A, and the eigenvectors of ATA correspond to the right singular vectors of A. If selected singular values and right singular vectors V' were computed, the related left singular vectors U' can be determined from the relation . The result is then processed with the modified Gram-Schmidt algorithm to assure that U' is orthogonal.

Similarly, since

the non-negative square roots of the eigenvalues of ATA correspond to the singular values of A, and the eigenvectors of AATcorrespond to the left singular vectors of A. For selected singular values and left singular vectors U', the related right singular vectors V' are computed from the relation and processed with the Gram-Schmidt algorithm to assure V' is orthogonal.

In the expanded approach, selectable by setting matrix_type = IMSL_ARPACK_EXPANDED, function imsl_d_arpack_svd iteratively computes selected eigenvalues and eigenvectors of the symmetric matrix

by applying function imsl_d_arpack_symmetric to this matrix internally. H is of order m+n and has eigenvalues and |m-n| additional zeros. If σi, i=1,…,r, are the positive eigenvalues of H, then a set of related orthonormal eigenvectors is

and a set of orthonormal eigenvectors for eigenvalue zero is

Therefore, the positive singular values and associated singular vectors of matrix A can be directly derived from the eigenvalues and eigenvectors of the expanded matrix H. The situation for eigenvalue zero and related eigenvectors is more difficult: First, if an eigenvalue zero is computed, it cannot be distinguished if it’s a singular value of A also or one of the |m-n| additional eigenvalues of value zero only. Second, splitting up eigenvectors in the same way as in the case of positive eigenvalues does not necessarily result in orthonormal singular vectors. As an example, consider matrix

Matrix

has the only eigenvalue 0 and a potential eigenbasis is

If the first two vectors are computed as eigenvectors, then, after identifying the first two components with a left singular vector and the last two components with a right singular vector,

but {u1,u2} and {v1,v2} are no orthogonal bases for the space of left and right singular vectors of A1.

Therefore, the expanded approach returns all selected and computable non-negative eigenvalues of H, but only singular vectors of A related to positive singular values.

If imsl_d_arpack_svd has been used to compute nsv singular values and associated left singular vectors U' and right singular vectors V', then the matrix B = U'Σ'V'T is an approximation to the input matrix A. Here,  Σ' is an nsv ×nsv diagonal matrix with the computed singular values in decreasing order on the diagonal. Especially, if =diag(σ1,σ2,…,σk) and are the k largest non-zero singular values of A, then the matrix

is the best rank -k approximation to A in the spectral norm .

The value of optional argument IMSL_MATRIX_TYPE may impact the number of iterations required. Generally, one expects matrix_type = IMSL_ARPACK_NORMAL (the default) to result in the fewest iterations, and matrix_type = IMSL_ARPACK_EXPANDED to result in singular values and vectors with the greatest accuracy.

Copyright Notice for ARPACK:

Copyright (c) 1996-2008 Rice University. Developed by D.C. Sorensen, R.B. Lehoucq, C. Yang, and K. Maschhoff. All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

- Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

- Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer listed in this license in the documentation and/or other materials provided with the distribution.

- Neither the name of the copyright holders nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Examples

Example 1

This example computes the two largest singular values of a real 6×4 matrix.

#include <imsl.h>

static void fcn(int m, int n, double x[], int itask, double y[]);

int main()
{
    int m = 6, n = 4, nsv = 2;
    int n_acc;
    double *s = NULL;

    /* Compute singular values */
    s = imsl_d_arpack_svd(fcn, m, n, nsv,
        IMSL_NUM_ACCURATE_SINGVALS, &n_acc,
        0);

    /* Print singular values */
    imsl_d_write_matrix("Largest singular values", 1, n_acc, s, 0);

    if (s)
        imsl_free(s);
}

static void fcn(int m, int n, double x[], int itask, double y[]) {

    double a[] = {
        1.0, 2.0, 1.0, 4.0,
        3.0, 2.0, 1.0, 3.0,
        4.0, 3.0, 1.0, 4.0,
        2.0, 1.0, 3.0, 1.0,
        1.0, 5.0, 2.0, 2.0,
        1.0, 2.0, 2.0, 3.0 };

    switch (itask) {
    case IMSL_ARPACK_PREPARE:
        /* Nothing to prepare, but formally handle this case */
        break;
    case IMSL_ARPACK_A_X:
        imsl_d_mat_mul_rect("A*x",
            IMSL_A_MATRIX, m, n, a,
            IMSL_X_VECTOR, n, x,
            IMSL_RETURN_USER, y,
            0);
        break;
    case IMSL_ARPACK_TRANS_A_X:
        imsl_d_mat_mul_rect("trans(A)*x",
            IMSL_A_MATRIX, m, n, a,
            IMSL_X_VECTOR, m, x,
            IMSL_RETURN_USER, y,
            0);
        break;
    default:
        imsl_set_user_fcn_return_flag(1);
        break;
    }
}

Output

Largest singular values
          1            2
      11.49         3.27

Example 2

We define the m×n matix A ={aij} by the entries aij=i+j. This matrix has at least two non-zero singular values. For the pair of values (m,n)= (600, 800) and (m,n)= (800, 600) we compute a partial SVD for the related rectangular matrices, using both matrix_type = IMSL_ARPACK_NORMAL and matrix_type = IMSL_ARPACK_EXPANDED. In each of the four possible cases, the resulting product B = UΣVT, when rounded to the nearest integer, satisfies B=A.

#include <imsl.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

/* Macro that rounds a real to the nearest integer */
#define round_to_int(x) (x < 0.0 ? (int)((x)-.5) : (int)((x)+.5))

typedef struct {
    double *a_matrix;
} imsl_arpack_data;

static void fcn_w_data(int m, int n, double x[], int itask, double y[],
    void *data);

static void compute_matrix_diffs(int m, int n, int nsv, double a[],
    double s[], double u[], double v[], double ab[]);

int main()
{
    int m = 600, n = 800, nsv = 2;
    int n_acc, max_m_n, count;
    double *a = NULL, *s = NULL, *u = NULL, *v = NULL, *ab = NULL;
    char *fmt = "%20.3f";
    imsl_arpack_data usr_data;

    /*
     *  Allocate memory for matrix "a" whose SVD will be computed
     *  by arpack_svd. The matrix will be filled with values under
     *  the IMSL_ARPACK_PREPARE task in user-defined function
     *  fcn_w_data().
     */
    a = (double *)malloc(m * sizeof(double) * n);
    usr_data.a_matrix = a;

    /*
     * The roles of m and n are interchanged in the following loop.
     * Accordingly, choose sizes of U and V so large that both cases
     * can be handled by the same arrays.
     */
    max_m_n = (m < n) ? n : m;
    u = (double *)malloc(max_m_n * sizeof(double) * nsv);
    v = (double *)malloc(max_m_n * sizeof(double) * nsv);

    /*
     * Allocate memory for singular value array "s" and auxiliary
     * array "ab".
     */
    s = (double *)malloc(nsv * sizeof(double));
    ab = (double *)malloc(m * sizeof(double) * n);

    for (count = 1; count <= 2; count++) {

        if (count == 2) {
            /* Interchange m and n */
            m = 800;
            n = 600;
        }

        printf("Matrix size: m = %d, n = %d\n", m, n);

        /* SVD computed with normal symmetric matrix */
        /* Compute singular values */
        imsl_d_arpack_svd(NULL, m, n, nsv,
            IMSL_SINGVAL_LOCATION, IMSL_ARPACK_LARGEST_ALGEBRAIC,
            IMSL_MATRIX_TYPE, IMSL_ARPACK_NORMAL,
            IMSL_NUM_ACCURATE_SINGVALS, &n_acc,
            IMSL_U_USER, u,
            IMSL_V_USER, v,
            IMSL_FCN_W_DATA, fcn_w_data, &usr_data,
            IMSL_RETURN_USER, s,
            0);

        /* Print singular values */
        imsl_d_write_matrix("Largest singular values, normal method",
            1, n_acc, s,
            IMSL_WRITE_FORMAT, fmt,
            0);
        printf("\nNumber of singular values, s, and columns of U and V : "
            "%d\n", n_acc);
        compute_matrix_diffs(m, n, n_acc, a, s, u, v, ab);

        /* SVD computed with expanded symmetric matrix */
        /* Compute singular values */
        imsl_d_arpack_svd(NULL, m, n, nsv,
            IMSL_SINGVAL_LOCATION, IMSL_ARPACK_LARGEST_ALGEBRAIC,
            IMSL_MATRIX_TYPE, IMSL_ARPACK_EXPANDED,
            IMSL_NUM_ACCURATE_SINGVALS, &n_acc,
            IMSL_U_USER, u,
            IMSL_V_USER, v,
            IMSL_FCN_W_DATA, fcn_w_data, &usr_data,
            IMSL_RETURN_USER, s,
            0);

        /* Print singular values */
        imsl_d_write_matrix("Largest singular values, expanded method",
            1, n_acc, s,
            IMSL_WRITE_FORMAT, fmt,
            0);
        printf("\nNumber of singular values, s, and columns of U and V : "
            "%d\n", n_acc);
        compute_matrix_diffs(m, n, n_acc, a, s, u, v, ab);
        if (count == 1) {
            printf("*****************************************************"
                "********\n\n");
        }
    }

    if (s)
        free(s);
    if (u)
        free(u);
    if (v)
        free(v);
    if (a)
        free(a);
    if (ab)
        free(ab);
}

static void fcn_w_data(int m, int n, double x[], int itask, double y[],
    void *data) {

    int i, j;
    double *a = NULL;

    imsl_arpack_data *usr_data = (imsl_arpack_data *)data;
    a = usr_data->a_matrix;

    switch (itask) {
    case IMSL_ARPACK_PREPARE:
        /* Fill matrix with integer values */
        for (i = 1; i <= m; i++) {
            for (j = 1; j <= n; j++) {
                a[(i - 1) * n + j - 1] = i + j;
            }
        }
        break;
    case IMSL_ARPACK_A_X:
        imsl_d_mat_mul_rect("A*x",
            IMSL_A_MATRIX, m, n, a,
            IMSL_X_VECTOR, n, x,
            IMSL_RETURN_USER, y,
            0);
        break;
    case IMSL_ARPACK_TRANS_A_X:
        imsl_d_mat_mul_rect("trans(A)*x",
            IMSL_A_MATRIX, m, n, a,
            IMSL_X_VECTOR, m, x,
            IMSL_RETURN_USER, y,
            0);
        break;
    default:
        imsl_set_user_fcn_return_flag(1);
        break;
    }
}

/*
 *  Calculates maximum absolute difference between integer matrix and
 *  rounded SVD approximation:
 *   max(abs(A - round(sum{u_i * sigma_i * trans(v_i) : i = 1,..., nsv})))
 */
static void compute_matrix_diffs(int m, int n, int nsv, double a[],
    double s[], double u[], double v[], double ab[]) {
    int i, j;
    double temp, max_diff;

    /* Compute W := U * S */
    for (i = 0; i < nsv; i++) {
        temp = s[i];
        for (j = 0; j < m; j++) {
            u[j * nsv + i] *= temp;
        }
    }

    /* Compute W * trans(V) = U * S * trans(V) */
    imsl_d_mat_mul_rect("A*trans(B)",
        IMSL_A_MATRIX, m, nsv, u,
        IMSL_B_MATRIX, n, nsv, v,
        IMSL_RETURN_USER, ab,
        0);

    /* Compute max(abs(A - round(U * S * trans(V)))) */
    max_diff = 0.0;
    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++) {
            ab[i * n + j] = a[i * n + j] - round_to_int(ab[i * n + j]);
            if (fabs(ab[i * n + j]) > max_diff) {
                max_diff = fabs(ab[i * n + j]);
            }
        }
    }

    if (max_diff > 0.0) {
        printf("\nmax_diff = %lf\n", max_diff);
        printf("Rounded approximation not identical with original "
            "integer matrix.\n\n");
    }
    else {
        printf("\nmax_diff = %lf\n", max_diff);
        printf("Rounded approximation identical with original "
            "integer matrix.\n\n");
    }
}

Output

Matrix size: m = 600, n = 800
 
  Largest singular values, normal method
                   1                     2
          523955.723             36644.238

Number of singular values, s, and columns of U and V : 2

max_diff = 0.000000
Rounded approximation identical with original integer matrix.

 
 Largest singular values, expanded method
                   1                     2
          523955.723             36644.238

Number of singular values, s, and columns of U and V : 2

max_diff = 0.000000
Rounded approximation identical with original integer matrix.

*************************************************************

Matrix size: m = 800, n = 600
 
  Largest singular values, normal method
                   1                     2
          523955.723             36644.238

Number of singular values, s, and columns of U and V : 2

max_diff = 0.000000
Rounded approximation identical with original integer matrix.

 
 Largest singular values, expanded method
                   1                     2
          523955.723             36644.238

Number of singular values, s, and columns of U and V : 2

max_diff = 0.000000
Rounded approximation identical with original integer matrix.

Warning Errors

For further warning errors, see function imsl_d_arpack_symmetric.

Fatal Errors

For further fatal errors, see function imsl_d_arpack_symmetric.