Linear algebra and embeddings

Settings and run statistics of the iterative eigensolvers arpack_rssolve and arpack_rnsolve.

typedef struct og_arpack_options_t {
    /* input */
    char bmat[1];        /* 'I': standard problem (the only kind supported) */
    int n;               /* order of the matrix */
    char which[2];       /* "LA", "SA", "LM", "SM", "BE" (symmetric);
        "LM", "SM", "LR", "SR", "LI", "SI" (general) */
    int nev;             /* number of eigenpairs wanted */
    og_real_t tol;       /* relative accuracy of the Ritz values; <= 0: machine eps */
    int ncv;             /* size of the search space; 0: chosen automatically */
    int ldv;             /* unused, keep 0 */
    int ishift;          /* unused, keep 1 */
    int mxiter;          /* maximum number of restarts */
    int nb;              /* unused, keep 1 */
    int mode;            /* 1: OP = A; 3: OP = inv(A - sigma I) (shift-invert) */
    int start;           /* 0: random start vector; 1: the caller's */
    int lworkl;          /* unused */
    og_real_t sigma;     /* shift for mode 3 */
    og_real_t sigmai;    /* imaginary part of the shift (general solver, unused) */
    /* output */
    int info;            /* 0: converged; 1: mxiter reached */
    int ierr;            /* always 0 */
    int noiter;          /* restarts taken */
    int nconv;           /* converged Ritz values */
    int numop;           /* products OP * x */
    int numopb;          /* unused, 0 */
    int numreo;          /* re-orthogonalisation passes */
    /* internal */
    int iparam[11];
    int ipntr[14];
} og_arpack_options_t;

The field names, types and order are those of the reference ARPACK interface, so existing code ports by renaming. Initialise it with arpack_options_init, then set at least n, nev and which. The solvers fill the output fields on return. Fields marked unused are kept for compatibility and ignored.

Parameters

• bmat 'I', a standard problem: the only kind supported (anything else is OG_UNIMPLEMENTED).
• n Order of the matrix.
• which Two characters choosing the wanted eigenvalues: LA, SA, LM, SM or BE for arpack_rssolve; LM, SM, LR, SR, LI or SI for arpack_rnsolve.
• nev Number of eigenvalues wanted.
• tol Relative accuracy of the Ritz values; 0 or less means machine precision.
• ncv Size of the Krylov search space; 0 chooses max(20, 2 nev + 1). It is raised to at least nev + 1 (nev + 2 for arpack_rnsolve) and cut to at most n.
• ldv Unused; keep 0.
• ishift Unused; keep 1.
• mxiter Maximum number of restarts (default 3000).
• nb Unused; keep 1.
• mode 1: the callback applies A. 3: shift-invert, the callback applies inv(A - sigma I).
• start 0: random start vector. 1: the caller's, from the first column of vectors (when it has n rows) or else storage->resid.
• lworkl Unused.
• sigma The shift for mode 3.
• sigmai Imaginary part of the shift; must be 0 (arpack_rnsolve refuses complex shifts).
• info Output: 0 when converged, 1 when mxiter was reached.
• ierr Output: always 0.
• noiter Output: restarts taken.
• nconv Output: number of converged Ritz values.
• numop Output: number of callback calls.
• numopb Output: unused, 0.
• numreo Output: re-orthogonalisation passes.
• iparam, ipntr Internal.

Checked

Same field names, types and order as the reference interface; the output fields info, noiter, nconv, numop and numreo are filled and tested.

C name

og_arpack_options_t

In notebooks

Og.ArpackOptions

igraph 1.0

igraph_arpack_options_t

See also

arpack_options_init, arpack_rssolve, arpack_rnsolve

Work memory for arpack_rssolve and arpack_rnsolve that can be reused across calls.

typedef struct og_arpack_storage_t {
    int maxn, maxncv, maxldv;
    og_real_t *v;
    og_real_t *workl;
    og_real_t *workd;
    og_real_t *d;
    og_real_t *resid;
    og_real_t *ax;
    int *select;
    og_real_t *di;
    og_real_t *workev;
} og_arpack_storage_t;

Made by arpack_storage_init for one kind of solver (symmetric or general) and freed by arpack_storage_destroy. maxn, maxncv and maxldv record the limits it was made for; the arrays are private workspaces of Open Graph's solvers and their layout is not ARPACK's. After a solve, resid holds the last residual direction, which a following call can use as its start vector (options->start = 1). A solver given NULL storage allocates its own.

Parameters

• maxn Largest matrix order the storage serves (read only).
• maxncv Largest search-space size it serves (read only).
• maxldv Recorded for compatibility (read only).
• v, workl, workd, d, resid, ax, select, di, workev Private work arrays; di and workev exist only for the general solver.

Checked

Same public fields as the reference interface; reuse across calls and refusal by the other solver kind are tested.

C name

og_arpack_storage_t

In notebooks

Og.ArpackStorage

igraph 1.0

igraph_arpack_storage_t

See also

arpack_storage_init, arpack_storage_destroy

The callback type through which the iterative eigensolvers multiply by the matrix.

typedef og_error_t og_arpack_function_t(
    og_real_t *to,
    const og_real_t *from,
    int n,
    void *extra
);

The solver calls fun(to, from, n, extra) with two arrays of length n; the callback writes OP * from into to, where OP is A in mode 1 and inv(A - sigma I) in mode 3. extra is the pointer passed to the solver. The matrix never has to be stored: any routine that can apply it will do. A return value other than OG_SUCCESS stops the solver, which then fails with that code.

Parameters

• to Output array of length n.
• from Input array of length n; do not modify it.
• n Length of the arrays.
• extra The caller's data, as passed to the solver.

Returns

OG_SUCCESS, or an error code that aborts the solve.

Checked

A non-success return stops the solver with that code (tested with OG_FAILURE).

C name

og_arpack_function_t

In notebooks

Og.ArpackFunction

igraph 1.0

igraph_arpack_function_t

See also

arpack_rssolve, arpack_rnsolve

Places the vertices in n dimensions using the leading eigenvectors (undirected) or singular vectors (directed) of the adjacency matrix.

og_error_t og_adjacency_spectral_embedding(
    const og_t *graph,
    og_int_t n,
    const og_vector_t *weights,
    og_eigen_which_position_t which,
    og_bool_t scaled,
    og_matrix_t *X,
    og_matrix_t *Y,
    og_vector_t *D,
    const og_vector_t *cvec,
    og_arpack_options_t *options
);

For an undirected graph the function takes n eigenpairs of A + diag(cvec), where A is the weighted adjacency matrix: parallel edges add their weights and a self-loop adds its weight once to the diagonal. which picks them: OG_EIGEN_LM the largest in magnitude, OG_EIGEN_LA the largest, OG_EIGEN_SA the smallest. X receives the eigenvectors as columns and D the eigenvalues; with scaled true each column is multiplied by sqrt(|eigenvalue|), the X = U |D|^1/2 of the adjacency spectral embedding (Sussman, Tang, Fishkind and Priebe 2012). Y, when given, receives a copy of X.

For a directed graph the function takes n singular triplets of A + diag(cvec) (OG_EIGEN_LM means the largest, as OG_EIGEN_LA): X from the left singular vectors, Y from the right ones, D the singular values; scaled multiplies both by sqrt(D). Right vectors of zero singular values are taken from the null space.

The eigenproblem is solved by arpack_rssolve. Each column's sign follows Open Graph's rule: its entry of largest magnitude is positive (for directed graphs Y follows X). When options is given its tol, ncv and mxiter are used and its output fields report the run; n, nev, which, mode and start are set by the function.

Parameters

• graph The graph.
• n Number of dimensions, 1 to |V|.
• weights One weight per edge, or NULL for 1. Negative weights are allowed.
• which OG_EIGEN_LM, OG_EIGEN_LA or OG_EIGEN_SA (enum og_eigen_which_position_t).
• scaled Whether to multiply the vectors by the square roots of the absolute values.
• X Initialised matrix; receives |V| x n coordinates.
• Y Initialised matrix or NULL. Directed graphs: the coordinates from the right singular vectors; undirected: a copy of X.
• D Initialised vector or NULL; receives the n eigenvalues (undirected) or singular values (directed).
• cvec Values added to the diagonal, one per vertex, or NULL for none.
• options Solver options or NULL for the defaults; on return its output fields describe the run.

Returns

OG_SUCCESS, or one of the error codes below.

Errors

• OG_EINVAL weights does not have one value per edge, n is outside 1..|V|, or which is not OG_EIGEN_LM, OG_EIGEN_LA or OG_EIGEN_SA.
• OG_EINVAL cvec does not have one value per vertex.
• OG_EARPACK The eigensolver did not converge within options->mxiter restarts.
• OG_ENOMEM Not enough memory.

Complexity

Each solver step multiplies once by the matrix, O(|V| + |E|) (twice for a directed graph), and orthogonalises, O(|V| x ncv); a restart adds O(|V| x ncv^2 + ncv^3). The number of restarts depends on how well the wanted values are separated from the rest.

Checked

Checked against igraph's own embedding code through R igraph 2.3.3 (igraph C 0.10) on 251 cases: values to 1e-8 of the largest (1e-6 directed), columns up to sign; also against python-igraph 1.0's eigen_adjacency, hub and authority scores and numpy. Column signs differ by design (igraph's follow its random start). An edgeless graph gives zeros where R igraph gives no values. Not tested against igraph 1.0 itself, which no binding exposes.

SIMD

Through the eigenvalue solvers (lapack_dsyevr, arpack_rssolve, arpack_rnsolve): their inner loops run on the SIMD kernels with a fixed order, the same bits on every CPU path.

C name

og_adjacency_spectral_embedding

In notebooks

Og.AdjacencySpectralEmbedding

igraph 1.0

igraph_adjacency_spectral_embedding

See also

laplacian_spectral_embedding, dim_select, arpack_rssolve

og_matrix_t X;
og_vector_t D;
og_matrix_init(&X, 0, 0);
og_vector_init(&D, 0);
/* two dimensions, largest eigenvalues in magnitude, scaled by sqrt(|D|) */
if (og_adjacency_spectral_embedding(&g, 2, NULL, OG_EIGEN_LM, true, &X, NULL, &D,
                                    NULL, NULL) == OG_SUCCESS) {
    /* row v of X holds the coordinates of vertex v */
}
og_vector_destroy(&D);
og_matrix_destroy(&X);

Places the vertices in n dimensions using the leading eigenvectors (undirected) or singular vectors (directed) of a Laplacian of the graph.

og_error_t og_laplacian_spectral_embedding(
    const og_t *graph,
    og_int_t n,
    const og_vector_t *weights,
    og_eigen_which_position_t which,
    og_laplacian_spectral_embedding_type_t type,
    og_bool_t scaled,
    og_matrix_t *X,
    og_matrix_t *Y,
    og_vector_t *D,
    og_arpack_options_t *options
);

type chooses the matrix. A is the adjacency matrix as in adjacency_spectral_embedding (a self-loop once on the diagonal) and D the diagonal matrix of strengths, the weighted degrees (an undirected self-loop counts twice): OG_EMBEDDING_D_A gives D - A, OG_EMBEDDING_DAD gives D^-1/2 A D^-1/2 and OG_EMBEDDING_I_DAD gives I - D^-1/2 A D^-1/2. A vertex of strength 0 gets 0 for its inverse square root. For directed graphs D - A uses the out-strengths, and DAD is O^-1/2 A P^-1/2 with O and P the out- and in-strengths.

Otherwise it works as adjacency_spectral_embedding: which (LM, LA or SA) selects eigenpairs (undirected) or singular triplets (directed), X, Y and D receive vectors and values, scaled multiplies by the square roots of the values, and each column's largest entry is made positive. The solver is arpack_rssolve.

Parameters

• graph The graph.
• n Number of dimensions, 1 to |V|.
• weights One weight per edge, or NULL for 1.
• which OG_EIGEN_LM, OG_EIGEN_LA or OG_EIGEN_SA (enum og_eigen_which_position_t).
• type OG_EMBEDDING_D_A, OG_EMBEDDING_DAD or OG_EMBEDDING_I_DAD (enum og_laplacian_spectral_embedding_type_t).
• scaled Whether to multiply the vectors by the square roots of the absolute values.
• X Initialised matrix; receives |V| x n coordinates.
• Y Initialised matrix or NULL. Directed graphs: the coordinates from the right singular vectors; undirected: a copy of X.
• D Initialised vector or NULL; receives the n eigenvalues or singular values.
• options Solver options or NULL for the defaults; on return its output fields describe the run.

Returns

OG_SUCCESS, or one of the error codes below.

Errors

• OG_EINVAL weights does not have one value per edge, n is outside 1..|V|, or which is not OG_EIGEN_LM, OG_EIGEN_LA or OG_EIGEN_SA.
• OG_EINVAL type is not one of the three og_laplacian_spectral_embedding_type_t values, or a strength is negative with OG_EMBEDDING_DAD or OG_EMBEDDING_I_DAD (possible with negative weights).
• OG_EARPACK The eigensolver did not converge within options->mxiter restarts.
• OG_ENOMEM Not enough memory.

Complexity

Each solver step multiplies once by the matrix, O(|V| + |E|) (twice for a directed graph), and orthogonalises, O(|V| x ncv); a restart adds O(|V| x ncv^2 + ncv^3). The number of restarts depends on how well the wanted values are separated from the rest.

Checked

Against igraph's own code through R igraph 2.3.3 (igraph C 0.10) on 270 cases: values to 1e-8 (1e-6 directed), columns up to sign; spectra of python-igraph 1.0's Laplacian matrices on 2002 more. Documented differences: weighted DAD and I_DAD follow the standard definition (igraph 0.10's values are not the spectrum of D^-1/2 W D^-1/2), a zero strength gives 0 where igraph gives NaN, and directed D_A and I_DAD, which igraph 0.10 refuses, are not tested against igraph.

SIMD

Through the eigenvalue solvers (lapack_dsyevr, arpack_rssolve, arpack_rnsolve): their inner loops run on the SIMD kernels with a fixed order, the same bits on every CPU path.

C name

og_laplacian_spectral_embedding

In notebooks

Og.LaplacianSpectralEmbedding

igraph 1.0

igraph_laplacian_spectral_embedding

See also

adjacency_spectral_embedding, dim_select, arpack_rssolve

Chooses how many leading values of a decreasing sequence (a scree plot) stand apart from the rest, by Zhu and Ghodsi's profile likelihood.

og_error_t og_dim_select(const og_vector_t *sv, og_int_t *dim);

For every split q = 1..p of the p values in sv into the first q values and the rest, the two groups are modelled as normal samples with their own means and a common variance (one sample when q = p), and the profile log-likelihood is computed (Zhu and Ghodsi 2006). The variance is the pooled within-group sum of squares over p - 2 degrees of freedom (p - 1 for one group). The q with the largest likelihood is returned. A split whose variance is zero or not finite is no candidate; when no split qualifies the answer is 0; a single value gives 1. Exact ties go to the smallest q.

The usual input is the eigenvalues or singular values of an embedding sorted in decreasing order; the function does not sort.

Parameters

• sv The values, normally in decreasing order; at least one.
• dim Receives the chosen number, 0 to p.

Returns

OG_SUCCESS, or one of the error codes below.

Errors

• OG_EINVAL sv is empty.
• OG_ENOMEM Not enough memory.

Complexity

O(p): one forward and one backward pass.

Checked

Equal to igraph's dim_select through R igraph 2.3.3 on 400 of 415 vectors. The other 15 have a split with exactly zero within-group variance, where igraph's answer follows its own rounding; on exact likelihood ties (about 15% of tie vectors) igraph sometimes takes the largest q, Open Graph always the smallest.

C name

og_dim_select

In notebooks

Og.DimSelect

igraph 1.0

igraph_dim_select

See also

adjacency_spectral_embedding, laplacian_spectral_embedding

Computes the dot product of two real vectors of equal length.

og_error_t og_blas_ddot(
    const og_vector_t *v1,
    const og_vector_t *v2,
    og_real_t *res
);

*res receives the sum of v1[i] * v2[i]. The arithmetic runs on Open Graph's SIMD kernels (AVX-512, AVX2, NEON or plain C, chosen at run time), which use one fixed order of operations, so the result is the same bit for bit on every CPU. The fixed order (16 partial sums joined by a fixed tree) can differ in the last bits from a plain left-to-right loop.

Parameters

• v1 The first vector.
• v2 The second vector, of the same length.
• res Receives the dot product.

Returns

OG_SUCCESS, or one of the error codes below.

Errors

• OG_EINVAL The vectors have different lengths.

Complexity

O(n).

Checked

200 random vector pairs (lengths 0 to 40, scales 1e-3 to 1e3) against numpy, to rounding.

SIMD

Runs on the SIMD dot-product kernel; the sum is taken in one fixed order on every path (AVX-512, AVX2, NEON, plain C), so every path gives the same bits.

C name

og_blas_ddot

In notebooks

Og.BlasDdot

igraph 1.0

igraph_blas_ddot

See also

blas_dnrm2, blas_dgemv

Returns the Euclidean norm of a real vector, without needless overflow or underflow.

og_real_t og_blas_dnrm2(const og_vector_t *v);

The sum of squares is first formed in one pass with the SIMD dot kernel; when it lies between 1e-280 and 1e280 its square root is the answer. Otherwise a second pass with a running scale factor computes the norm, so vectors with entries near 1e200 or 1e-200 still get the right value. A NaN entry gives NaN and an infinite one infinity; an empty vector has norm 0. The result is the same on every CPU path.

Parameters

• v The vector.

Returns

The norm. The function cannot fail.

Complexity

O(n).

Checked

200 vectors, including entries near 1e-200 and 1e200, against math.hypot.

SIMD

The sum of squares runs on the SIMD dot-product kernel, in one fixed order on every path, so every path gives the same bits.

C name

og_blas_dnrm2

In notebooks

Og.BlasDnrm2

igraph 1.0

igraph_blas_dnrm2

See also

blas_ddot

Computes y = alpha A x + beta y, or the same with A transposed.

og_error_t og_blas_dgemv(
    og_bool_t transpose,
    og_real_t alpha,
    const og_matrix_t *a,
    const og_vector_t *x,
    og_real_t beta,
    og_vector_t *y
);

A is m x n. With transpose false x must have n entries and y m; with transpose true the other way round. When beta is 0 the old contents of y are never read (NaN there does no harm); when alpha is 0, A and x are not used. As in the reference BLAS, a column whose multiplier alpha x[j] is zero is skipped, so infinities or NaN in it do not reach y. The arithmetic runs on Open Graph's SIMD kernels (AVX-512, AVX2, NEON or plain C, chosen at run time), which use one fixed order of operations, so the result is the same bit for bit on every CPU.

Parameters

• transpose Whether to use the transpose of A.
• alpha Factor of the product.
• a The matrix A.
• x The input vector.
• beta Factor of the old y.
• y Input and output vector.

Returns

OG_SUCCESS, or one of the error codes below.

Errors

• OG_EINVAL x or y has the wrong length.
• OG_EOVERFLOW A dimension of a matrix is above INT_MAX (kept from the reference BLAS interface, which passes 32-bit sizes).

Complexity

O(m n).

Checked

200 cases (0 to 7 rows and columns, both transposes, alpha and beta including 0) against numpy.

SIMD

The row or column products run on the SIMD kernels with a fixed summation order; the same bits on every CPU path.

C name

og_blas_dgemv

In notebooks

Og.BlasDgemv

igraph 1.0

igraph_blas_dgemv

See also

blas_dgemv_array, blas_dgemm

Computes C = alpha op(A) op(B) + beta C, where op is the matrix itself or its transpose.

og_error_t og_blas_dgemm(
    og_bool_t transpose_a,
    og_bool_t transpose_b,
    og_real_t alpha,
    const og_matrix_t *a,
    const og_matrix_t *b,
    og_real_t beta,
    og_matrix_t *c
);

op(A) must be m x k and op(B) k x n. When beta is 0, c is resized to m x n and its old contents are ignored; otherwise c must already be m x n. c may be the same matrix as a or b: the product is then formed in a temporary copy. A zero multiplier skips its column as in the reference BLAS. Rows are processed in blocks of 256 so a tall A stays in cache. The arithmetic runs on Open Graph's SIMD kernels (AVX-512, AVX2, NEON or plain C, chosen at run time), which use one fixed order of operations, so the result is the same bit for bit on every CPU.

Parameters

• transpose_a Whether to use the transpose of A.
• transpose_b Whether to use the transpose of B.
• alpha Factor of the product.
• a The matrix A.
• b The matrix B.
• beta Factor of the old C.
• c Input and output matrix.

Returns

OG_SUCCESS, or one of the error codes below.

Errors

• OG_EINVAL The inner dimensions differ, or beta is not 0 and c is not m x n.
• OG_EOVERFLOW A dimension of a matrix is above INT_MAX (kept from the reference BLAS interface, which passes 32-bit sizes).
• OG_ENOMEM Not enough memory.

Complexity

O(m n k).

Checked

200 cases (both transposes, empty dimensions, alpha and beta including 0) against numpy.

SIMD

The inner products run on the SIMD kernels with a fixed summation order; the same bits on every CPU path.

C name

og_blas_dgemm

In notebooks

Og.BlasDgemm

igraph 1.0

igraph_blas_dgemm

See also

blas_dgemv

The product of blas_dgemv with plain arrays for x and y.

og_error_t og_blas_dgemv_array(
    og_bool_t transpose,
    og_real_t alpha,
    const og_matrix_t *a,
    const og_real_t *x,
    og_real_t beta,
    og_real_t *y
);

Computes y = alpha op(A) x + beta y exactly as blas_dgemv does, but x and y are raw arrays whose lengths are not checked: x must hold at least ncol(A) values (nrow(A) when transposed) and y at least nrow(A) (ncol(A) when transposed).

Parameters

• transpose Whether to use the transpose of A.
• alpha Factor of the product.
• a The matrix A.
• x The input array.
• beta Factor of the old y.
• y Input and output array.

Returns

OG_SUCCESS, or one of the error codes below.

Errors

• OG_EOVERFLOW A dimension of a matrix is above INT_MAX (kept from the reference BLAS interface, which passes 32-bit sizes).

Complexity

O(m n).

Checked

The blas_dgemv fixture run through the array version.

SIMD

As blas_dgemv(): the SIMD kernels with a fixed summation order, the same bits on every CPU path.

C name

og_blas_dgemv_array

In notebooks

Og.BlasDgemvArray

igraph 1.0

igraph_blas_dgemv_array

See also

blas_dgemv

Factors a general m x n matrix as P A = L U by Gaussian elimination with partial pivoting, in place.

og_error_t og_lapack_dgetrf(
    og_matrix_t *a,
    og_vector_int_t *ipiv,
    int *info
);

On return a holds U on and above the diagonal and the multipliers of L (unit lower triangular; its diagonal is not stored) below it. ipiv receives min(m, n) pivot rows, 1-based as in LAPACK: row k was swapped with row ipiv[k]. *info is 0, or k > 0 when U(k,k) is exactly zero (the first such k); the factorisation is still completed, but U is singular and must not be used to solve.

The multipliers are formed by division rather than by multiplication with a reciprocal, so a row equal to the pivot row cancels to an exact zero on every platform.

Parameters

• a The matrix; overwritten by the factors.
• ipiv Receives the 1-based pivot rows; may be NULL.
• info Receives 0 or the 1-based index of the first zero pivot; may be NULL.

Returns

OG_SUCCESS, or one of the error codes below.

Errors

• OG_ENOMEM Not enough memory.

Complexity

O(m n min(m, n)); about 2n^3/3 flops for a square matrix.

Checked

Pivots and info identical to scipy's LAPACK on 200 matrices (square, rectangular, singular), factors to 1e-10. A pivot that cancels to zero only through rounding can be exactly 0 in OpenBLAS (fused multiply-add) and 1e-16 here, so info may differ there.

C name

og_lapack_dgetrf

In notebooks

Og.LapackDgetrf

igraph 1.0

igraph_lapack_dgetrf

See also

lapack_dgetrs, lapack_dgesv

Solves A X = B or A^T X = B with the LU factors from lapack_dgetrf.

og_error_t og_lapack_dgetrs(
    og_bool_t transpose,
    const og_matrix_t *a,
    const og_vector_int_t *ipiv,
    og_matrix_t *b
);

a and ipiv must be the output of lapack_dgetrf for a square matrix of order n. b holds the right-hand sides as its n x nrhs columns and is overwritten by the solutions; transpose true solves the transposed system. Zero pivots are not checked (dividing by one gives infinities or NaN), so check info from the factorisation first.

Parameters

• transpose Whether to solve with A^T.
• a The LU factors.
• ipiv The 1-based pivots.
• b The right-hand sides; overwritten by the solutions.

Returns

OG_SUCCESS, or one of the error codes below.

Errors

• OG_EINVAL a is not square, b does not have n rows, ipiv does not have n entries, or a pivot lies outside 1..n.
• OG_ENOMEM Not enough memory.

Complexity

O(n^2 x nrhs).

Checked

200 systems (0 to 9 unknowns, 0 to 3 right-hand sides, some singular), plain and transposed, against numpy.linalg.solve.

C name

og_lapack_dgetrs

In notebooks

Og.LapackDgetrs

igraph 1.0

igraph_lapack_dgetrs

See also

lapack_dgetrf, lapack_dgesv

Solves the square linear system A X = B by LU factorisation with partial pivoting.

og_error_t og_lapack_dgesv(
    og_matrix_t *a,
    og_vector_int_t *ipiv,
    og_matrix_t *b,
    int *info
);

Factors a in place as lapack_dgetrf does and, when no pivot is exactly zero, overwrites b with the solutions. ipiv receives the 1-based pivots. When a pivot is zero, *info is its 1-based index, a warning is given and b is left unchanged; the call still returns OG_SUCCESS, so check info.

Parameters

• a The square matrix; overwritten by its LU factors.
• ipiv Receives the 1-based pivots; may be NULL.
• b The right-hand sides (n x nrhs); overwritten by the solutions.
• info Receives 0 or the index of the first zero pivot; may be NULL.

Returns

OG_SUCCESS, or one of the error codes below.

Errors

• OG_EINVAL a is not square, or b does not have n rows.
• OG_ENOMEM Not enough memory.

Complexity

O(n^3 + n^2 x nrhs).

Checked

The same 200 systems against numpy.linalg.solve; info exactly as scipy's LAPACK dgetrf.

C name

og_lapack_dgesv

In notebooks

Og.LapackDgesv

igraph 1.0

igraph_lapack_dgesv

See also

lapack_dgetrf, lapack_dgetrs

Eigen-decomposition of a real symmetric matrix: every eigenvalue or a chosen subset, with their eigenvectors when asked for.

og_error_t og_lapack_dsyevr(
    const og_matrix_t *A,
    og_lapack_dsyev_which_t which,
    og_real_t vl,
    og_real_t vu,
    int vestimate,
    int il,
    int iu,
    og_real_t abstol,
    og_vector_t *values,
    og_matrix_t *vectors,
    og_vector_int_t *support
);

Only the upper triangle of A is read. The matrix is reduced to tridiagonal form by Householder reflections and diagonalised by the implicit symmetric QR iteration with Wilkinson shifts (not LAPACK's MRRR algorithm: the values agree to rounding, the vectors up to sign). which selects the output: OG_LAPACK_DSYEV_ALL every eigenvalue, OG_LAPACK_DSYEV_INTERVAL those in the half-open interval (vl, vu], OG_LAPACK_DSYEV_SELECT those with 1-based positions il..iu in ascending order. values receives them in ascending order and vectors the matching orthonormal eigenvectors as columns.

support, when given, receives two 1-based row indices per eigenvector: its first and last non-zero row. vestimate and abstol are accepted for compatibility and not used; the outputs are always sized exactly. The arithmetic runs on Open Graph's SIMD kernels (AVX-512, AVX2, NEON or plain C, chosen at run time), which use one fixed order of operations, so the result is the same bit for bit on every CPU.

Parameters

• A The symmetric matrix; only its upper triangle is read.
• which OG_LAPACK_DSYEV_ALL, OG_LAPACK_DSYEV_INTERVAL or OG_LAPACK_DSYEV_SELECT (enum og_lapack_dsyev_which_t).
• vl INTERVAL: lower bound, excluded.
• vu INTERVAL: upper bound, included.
• vestimate Accepted for compatibility; not used.
• il SELECT: 1-based position of the first eigenvalue wanted.
• iu SELECT: 1-based position of the last eigenvalue wanted.
• abstol Accepted for compatibility; not used.
• values Receives the eigenvalues in ascending order; may be NULL.
• vectors Receives the eigenvectors as columns; may be NULL.
• support Receives the first and last non-zero row of each eigenvector (1-based); may be NULL.

Returns

OG_SUCCESS, or one of the error codes below.

Errors

• OG_EINVAL A is not square, vl >= vu with INTERVAL, il..iu is not a range inside 1..n with SELECT, or which is unknown.
• OG_DIVERGED The QR iteration did not converge.
• OG_ENOMEM Not enough memory.

Complexity

O(n^3).

Checked

220 symmetric matrices up to 64 x 64 (dense, adjacency, Laplacian, diagonal, zero) with ALL, INTERVAL and SELECT against scipy's LAPACK dsyevr: values to 1e-10, vectors up to sign, orthonormality and residuals. support is each vector's true non-zero range, not LAPACK's isuppz.

SIMD

The inner loops run on the SIMD kernels with a fixed order; the same bits on every CPU path.

C name

og_lapack_dsyevr

In notebooks

Og.LapackDsyevr

igraph 1.0

igraph_lapack_dsyevr

See also

lapack_dgeev, arpack_rssolve

og_vector_t values;
og_matrix_t vectors;
og_vector_init(&values, 0);
og_matrix_init(&vectors, 0, 0);
/* the two smallest eigenpairs of the symmetric matrix A */
og_lapack_dsyevr(&A, OG_LAPACK_DSYEV_SELECT, 0, 0, 0, 1, 2, 0, &values, &vectors, NULL);
og_matrix_destroy(&vectors);
og_vector_destroy(&values);

Computes the eigenvalues and, optionally, the left and right eigenvectors of a general real square matrix.

og_error_t og_lapack_dgeev(
    const og_matrix_t *A,
    og_vector_t *valuesreal,
    og_vector_t *valuesimag,
    og_matrix_t *vectorsleft,
    og_matrix_t *vectorsright,
    int *info
);

The matrix is balanced (permutation and scaling), reduced to upper Hessenberg form and brought to real Schur form by the Francis double-shift QR iteration; the eigenvectors come from back substitution and are transformed back. The real and imaginary parts go to valuesreal and valuesimag; a complex conjugate pair takes two adjacent places, the member with positive imaginary part first. In vectorsright and vectorsleft a real eigenvalue has one column and a conjugate pair two: the real and imaginary parts of the eigenvector of the first member (the second member's is its conjugate). Every eigenvector has unit 2-norm with its largest component real and positive.

*info on entry decides what a failure of the QR iteration means: non-zero (or info NULL) makes it the error OG_DIVERGED; zero makes it a warning, with the converged eigenvalues returned and the eigenvectors set to NaN. On exit *info is 0, or the number of leading eigenvalues that were not found (positions info+1..n are valid). Only the eigenvector norms use the SIMD kernels; the iteration itself is plain C.

Parameters

• A The square matrix; not modified.
• valuesreal Receives the real parts of the eigenvalues; may be NULL.
• valuesimag Receives the imaginary parts; may be NULL.
• vectorsleft Receives the left eigenvectors (n x n); may be NULL.
• vectorsright Receives the right eigenvectors (n x n); may be NULL.
• info On entry: non-zero to treat non-convergence as an error. On exit: 0, or the count of eigenvalues not found.

Returns

OG_SUCCESS, or one of the error codes below.

Errors

• OG_EINVAL A is not square, or has a NaN or infinite entry.
• OG_DIVERGED The QR iteration did not converge and info asked for an error.
• OG_ENOMEM Not enough memory.

Complexity

O(n^3).

Checked

Same values and vectors as lapack_dgeevx with both balancings, which is checked against LAPACK; vectors compared up to phase (LAPACK leaves the sign of a real vector to its back substitution). Hard cases (permutation matrices, Jordan blocks, graded and random orthogonal matrices, 204 matrices with repeated eigenvalues) converge with residuals below 1e-13. The non-convergence path is not tested: no matrix was found that makes the iteration fail.

SIMD

The inner loops run on the SIMD kernels with a fixed order; the same bits on every CPU path.

C name

og_lapack_dgeev

In notebooks

Og.LapackDgeev

igraph 1.0

igraph_lapack_dgeev

See also

lapack_dgeevx, lapack_dsyevr, arpack_rnsolve

lapack_dgeev with a choice of balancing, the balancing details and reciprocal condition numbers.

og_error_t og_lapack_dgeevx(
    og_lapack_dgeevx_balance_t balance,
    const og_matrix_t *A,
    og_vector_t *valuesreal,
    og_vector_t *valuesimag,
    og_matrix_t *vectorsleft,
    og_matrix_t *vectorsright,
    int *ilo,
    int *ihi,
    og_vector_t *scale,
    og_real_t *abnrm,
    og_vector_t *rconde,
    og_vector_t *rcondv,
    int *info
);

balance chooses the balancing: OG_LAPACK_DGEEVX_BALANCE_NONE, _PERM (permute rows and columns to isolate eigenvalues), _SCALE (diagonal scaling by powers of two to even out row and column norms) or _BOTH, which is what lapack_dgeev uses. ilo and ihi receive the 1-based range of rows and columns the permutation left in place; scale receives, for each row outside ilo..ihi, the 1-based row it was exchanged with, and inside the range its scaling factor; abnrm receives the 1-norm (largest column sum of absolute values) of the balanced matrix.

rconde receives the reciprocal condition number of each eigenvalue and rcondv that of each right eigenvector, both of the balanced matrix as defined in the LAPACK Users' Guide; rcondv uses an estimate of sep by the Hager-Higham 1-norm estimator. Small values mean sensitive results. Values, vectors and info work as in lapack_dgeev.

Parameters

• balance OG_LAPACK_DGEEVX_BALANCE_NONE, _PERM, _SCALE or _BOTH (enum og_lapack_dgeevx_balance_t).
• A The square matrix; not modified.
• valuesreal Receives the real parts of the eigenvalues; may be NULL.
• valuesimag Receives the imaginary parts; may be NULL.
• vectorsleft Receives the left eigenvectors (n x n); may be NULL.
• vectorsright Receives the right eigenvectors (n x n); may be NULL.
• ilo Receives the first row of the balanced block (1-based); may be NULL.
• ihi Receives the last row of the balanced block (1-based); may be NULL.
• scale Receives the permutation and scaling details, one entry per row; may be NULL.
• abnrm Receives the 1-norm of the balanced matrix; may be NULL.
• rconde Receives the reciprocal condition numbers of the eigenvalues; may be NULL.
• rcondv Receives the reciprocal condition numbers of the right eigenvectors; may be NULL.
• info As in lapack_dgeev.

Returns

OG_SUCCESS, or one of the error codes below.

Errors

• OG_EINVAL balance is unknown, or A is not square or has a NaN or infinite entry.
• OG_DIVERGED The QR iteration did not converge and info asked for an error.
• OG_ENOMEM Not enough memory.

Complexity

O(n^3).

Checked

270 matrices x 4 balancings against LAPACK's DGEEVX called as a black box: ilo, ihi and scale identical bit for bit, abnrm to 1e-13, eigenvalues within their condition (in LAPACK's order in 264 of 270), vectors up to phase, rconde to 1e-7 where well conditioned; rcondv (a 1-norm estimate in both libraries) agrees to 1e-6 in 1446 of 1709 comparisons and within the estimate's bound everywhere.

SIMD

As lapack_dgeev(): the SIMD kernels with a fixed order, the same bits on every CPU path.

C name

og_lapack_dgeevx

In notebooks

Og.LapackDgeevx

igraph 1.0

igraph_lapack_dgeevx

See also

lapack_dgeev

Fills an arpack_options_t with the default settings.

void og_arpack_options_init(og_arpack_options_t *o);

Sets bmat 'I', nev 1, tol 0 (machine precision), ncv 0 (automatic), ishift 1, mxiter 3000, nb 1, mode 1, start 0 (random start vector) and every other field to 0. which is set to "XX", which the solvers refuse with OG_EINVAL, so the caller must choose it, as well as n and nev.

Parameters

• o The options to fill.

Returns

Nothing.

Complexity

O(1).

Checked

The documented defaults; an unset which is refused by the solvers.

C name

og_arpack_options_init

In notebooks

Og.ArpackOptionsInit

igraph 1.0

igraph_arpack_options_init

See also

arpack_options_t, arpack_rssolve, arpack_rnsolve

Allocates reusable work memory for the iterative eigensolvers.

og_error_t og_arpack_storage_init(
    og_arpack_storage_t *s,
    og_int_t maxn,
    og_int_t maxncv,
    og_int_t maxldv,
    og_bool_t symm
);

maxn is the largest matrix order and maxncv the largest search-space size the storage must serve; maxldv is only recorded. symm true makes storage for arpack_rssolve, false for arpack_rnsolve, which needs two more arrays; a solver refuses storage made for the other kind or too small for its n and ncv. The memory is about 2 maxn (maxncv + 1) + 6 maxncv^2 reals. On failure nothing stays allocated.

Parameters

• s The storage to initialise.
• maxn Largest matrix order.
• maxncv Largest search-space size.
• maxldv Recorded only.
• symm true for the symmetric solver, false for the general one.

Returns

OG_SUCCESS, or one of the error codes below.

Errors

• OG_EINVAL A size is negative or above INT_MAX.
• OG_EOVERFLOW The arrays would be too large.
• OG_ENOMEM Not enough memory.

Complexity

O(maxn x maxncv + maxncv^2) (the memory is zeroed).

Checked

Size checks, refusal of too small or wrong-kind storage, and the out-of-memory path are unit-tested.

C name

og_arpack_storage_init

In notebooks

Og.ArpackStorageInit

igraph 1.0

igraph_arpack_storage_init

See also

arpack_storage_t, arpack_storage_destroy

Frees the memory of an arpack_storage_t.

void og_arpack_storage_destroy(og_arpack_storage_t *s);

Frees every array and sets the pointers to NULL, so a second call on the same storage is harmless.

Parameters

• s The storage to free.

Returns

Nothing.

Complexity

O(1).

Checked

Safe on an already destroyed storage.

C name

og_arpack_storage_destroy

In notebooks

Og.ArpackStorageDestroy

igraph 1.0

igraph_arpack_storage_destroy

See also

arpack_storage_init

Finds a few eigenvalues and eigenvectors of a real symmetric matrix that is given only as a matrix-vector product.

og_error_t og_arpack_rssolve(
    og_arpack_function_t *fun,
    void *extra,
    og_arpack_options_t *options,
    og_arpack_storage_t *storage,
    og_vector_t *values,
    og_matrix_t *vectors
);

Open Graph's own solver behind the ARPACK interface; no ARPACK code is used. The method is thick-restart Lanczos (Wu and Simon 2000) with full re-orthogonalisation. The callback applies the matrix (mode 1) or inv(A - sigma I) (mode 3, shift-invert; the values are mapped back to sigma + 1/theta). options->which chooses the eigenvalues: LA largest, SA smallest, LM largest in magnitude, SM smallest in magnitude, BE from both ends. values receives nev of them in that order (LA descending, SA ascending, LM and SM by magnitude with the larger value first on ties, BE alternating largest, smallest, ...) and vectors the n x nev orthonormal eigenvectors.

A Krylov space holds only one copy of a repeated eigenvalue. After convergence the solver therefore restarts from a random direction orthogonal to the converged vectors (twice in a row when the problem has shown a repeated or zero eigenvalue) and accepts the set only when nothing better turns up; ARPACK can return an incomplete set there. Random start vectors come from the default random number generator.

Orthogonalisation and the basis updates run on the SIMD kernels with a fixed order of operations.

Parameters

• fun The callback that applies the matrix (arpack_function_t).
• extra Passed to fun unchanged.
• options Settings on entry (n, nev, which, ...); run statistics on return.
• storage Work memory from arpack_storage_init(..., true), or NULL to allocate it internally.
• values Receives the nev eigenvalues; may be NULL.
• vectors Receives the n x nev eigenvectors; may be NULL. With options->start = 1 its first column is the start vector when it has n rows.

Returns

OG_SUCCESS, one of the error codes below, or the error code the callback returned (the solver stops at once).

Errors

• OG_EINVAL n is negative, which is not LA, SA, LM, SM or BE, nev is outside 1..n, a start vector is requested but none is given, or storage is too small or made for the general solver.
• OG_UNIMPLEMENTED bmat is not 'I', or mode is neither 1 nor 3.
• OG_EARPACK Fewer than nev values converged within mxiter restarts (options->info is then 1), the start vector is zero or not finite, or the projected problem failed.
• OG_ENOMEM Not enough memory.

Complexity

Per restart: up to ncv callback calls, O(n x ncv^2) for re-orthogonalisation and O(ncv^3) for the projected problem. The number of restarts depends on the gap between wanted and unwanted eigenvalues; mxiter bounds it.

Checked

200 undirected random graphs (weights, loops, multi-edges) with all five codes and two ncv settings against numpy's dense eigenvalues: values to 1e-8 ||A||, residuals and orthonormality; 27 cases also against python-igraph's eigen_adjacency. In regular mode with SM on interior zero eigenvalues of high multiplicity scipy's ARPACK fails 3 fixture cases; Open Graph answers 2 and stops with OG_EARPACK on 1.

SIMD

The Lanczos solver's dot products and vector updates run on the SIMD kernels; reductions use one fixed order, so every CPU path gives the same bits.

C name

og_arpack_rssolve

In notebooks

Og.ArpackRssolve

igraph 1.0

igraph_arpack_rssolve

See also

arpack_rnsolve, arpack_options_t, arpack_function_t, lapack_dsyevr

/* to = A from, for a dense symmetric matrix passed in extra */
static og_error_t apply(og_real_t *to, const og_real_t *from, int n, void *extra) {
    const og_matrix_t *A = (const og_matrix_t *) extra;
    for (int i = 0; i < n; i++) {
        og_real_t s = 0.0;
        for (int j = 0; j < n; j++) {
            s += MATRIX(*A, i, j) * from[j];
        }
        to[i] = s;
    }
    return OG_SUCCESS;
}

/* ... the three largest eigenpairs ... */
og_arpack_options_t opt;
og_vector_t values;
og_matrix_t vectors;
og_arpack_options_init(&opt);
opt.n = (int) og_matrix_nrow(&A);
opt.nev = 3;
opt.which[0] = 'L';
opt.which[1] = 'A';
og_vector_init(&values, 0);
og_matrix_init(&vectors, 0, 0);
og_arpack_rssolve(apply, &A, &opt, NULL, &values, &vectors);

Finds a few eigenvalues and eigenvectors of a real non-symmetric matrix that is given only as a matrix-vector product.

og_error_t og_arpack_rnsolve(
    og_arpack_function_t *fun,
    void *extra,
    og_arpack_options_t *options,
    og_arpack_storage_t *storage,
    og_matrix_t *values,
    og_matrix_t *vectors
);

Open Graph's own solver behind the ARPACK interface; no ARPACK code is used. The method is Krylov-Schur restarted Arnoldi (Stewart 2001), equivalent to implicit restarting with exact shifts. options->which selects LM or SM (largest or smallest magnitude), LR or SR (largest or smallest real part), LI or SI (largest or smallest imaginary part). When the search space would reach n (small matrices), the matrix is built column by column through the callback and solved densely. Mode 3 (shift-invert with a real sigma) maps the values back to sigma + 1/theta.

values receives a k x 2 matrix (real part, imaginary part), where k is nev, or nev + 1 when the nev-th value is the first member of a conjugate pair, so pairs are never split; pairs are adjacent with the positive imaginary part first. vectors receives n x k columns: a real eigenvalue takes one column, a conjugate pair two (the real and imaginary parts of the eigenvector of the member with positive imaginary part). arpack_unpack_complex turns this into one complex vector per value. Random start vectors come from the default random number generator; orthogonalisation runs on the SIMD kernels.

Parameters

• fun The callback that applies the matrix (arpack_function_t).
• extra Passed to fun unchanged.
• options Settings on entry (n, nev, which, ...); run statistics on return.
• storage Work memory from arpack_storage_init(..., false), or NULL to allocate it internally.
• values Receives the k x 2 eigenvalues; may be NULL.
• vectors Receives the n x k eigenvector columns; may be NULL. With options->start = 1 its first column is the start vector when it has n rows.

Returns

OG_SUCCESS, one of the error codes below, or the error code the callback returned (the solver stops at once).

Errors

• OG_EINVAL n is negative, which is not LM, SM, LR, SR, LI or SI, nev is outside 1..n, a start vector is requested but none is given, or storage is too small or made for the symmetric solver.
• OG_UNIMPLEMENTED sigmai is not 0 in mode 3 (complex shifts).
• OG_UNIMPLEMENTED bmat is not 'I', or mode is neither 1 nor 3.
• OG_EARPACK Fewer than nev values converged within mxiter restarts (options->info is then 1), the start vector is zero or not finite, or the projected problem failed.
• OG_ENOMEM Not enough memory.

Complexity

Per restart: up to ncv callback calls, O(n x ncv^2) for orthogonalisation and O(ncv^3) for the Schur form; when ncv reaches n, n callback calls and a dense O(n^3) solve.

Checked

260 directed adjacency matrices (signed weights included) with all six codes against the full spectrum from LAPACK: the best values wherever scipy's ARPACK finds them (174 cases) and 36 of the 42 where it fails; 65 Perron roots equal python-igraph's eigenvector_centrality eigenvalue to 1e-8; 96 dense random matrices and shift-invert mode. Not matched: defective zero eigenvalues (nilpotent parts of DAG-like graphs) come back at their sensitivity, about 1e-6, where ARPACK, iterating longer, often lands on 0.

SIMD

The Arnoldi solver's dot products and vector updates run on the SIMD kernels; reductions use one fixed order, so every CPU path gives the same bits.

C name

og_arpack_rnsolve

In notebooks

Og.ArpackRnsolve

igraph 1.0

igraph_arpack_rnsolve

See also

arpack_unpack_complex, arpack_rssolve, lapack_dgeev

Rewrites arpack_rnsolve's output so that every eigenvalue has its own complex eigenvector.

og_error_t og_arpack_unpack_complex(
    og_matrix_t *vectors,
    og_matrix_t *values,
    og_int_t nev
);

values (k x 2) and vectors are replaced in place: values keeps its first nev rows, and vectors becomes n x 2 nev, where columns 2i and 2i+1 hold the real and imaginary parts of the eigenvector of value i. A real value gets a zero imaginary column; the second member of a conjugate pair gets the conjugate of the first member's vector.

Parameters

• vectors arpack_rnsolve's eigenvectors; replaced by n x 2 nev columns.
• values arpack_rnsolve's k x 2 eigenvalues; cut to nev rows.
• nev Number of eigenvalues to keep, 0 to k.

Returns

OG_SUCCESS, or one of the error codes below.

Errors

• OG_EINVAL nev is negative or above the number of values, values does not have two columns, or vectors has too few columns for the values.
• OG_ENOMEM Not enough memory.

Complexity

O(n x nev).

Checked

Hand-built inputs (real values, pairs, a lone member with negative imaginary part, fewer values kept) and arpack_rnsolve's own output; python-igraph has no counterpart.

C name

og_arpack_unpack_complex

In notebooks

Og.ArpackUnpackComplex

igraph 1.0

igraph_arpack_unpack_complex

See also

arpack_rnsolve