Ritz values (Rayleigh quotients)

Status: current · signatures verified against include/batchlas/blas/extensions.hh (lines 1621-1672) and src/ops/ritz_values/ritz_values.cc on 2026-09-30.

ritz_values returns the Rayleigh quotient of every trial-vector column, for dense or CSR operators and batched inputs.

What it computes

For square \(A\) (dense or CSR) and a dense block \(V = [v_1, \dots, v_k]\), per column and batch item:

\[ \theta_j \;=\; \operatorname{Re}\,\frac{v_j^H A\, v_j}{v_j^H v_j}, \qquad j = 1, \dots, k . \]

• \(v_j^H = v_j^T\) for real types. The quotient is formed in the full scalar type and only the real part is stored, so the output is float for std::complex<float> and double for std::complex<double>.
• For Hermitian \(A\) and \(v_j\) near an eigenvector, \(\theta_j\) is second-order accurate. stedc_tests uses it to recompute eigenvalues from returned eigenvectors; Python exposes it as batchlas.ritz_values(a, vectors).
• Columns are independent: the result is \(\operatorname{diag}(V^H A V)\) divided by \(\operatorname{diag}(V^H V)\), not the eigenvalues of \(V^H A V\). Columns need not be normalised or orthogonal.

Note: For a Rayleigh-Ritz projection, form \(V^H A V\) with gemm and call syev on it.

How it is computed

1. \(W = A V\) into workspace: gemm for dense \(A\), spmm for CSR (its workspace is carved from the same buffer).
2. One kernel, one work-group per (batch item, column) pair, grid-strided. Each work-group reduces \(v_j^H w_j\) and \(v_j^H v_j\) through local memory; work-item 0 writes \(\operatorname{Re}(\text{num}/\text{den})\).

Calling it

Allocating overload

Allocates the output and workspace itself and blocks (ctx.wait() before and after). Keep it out of loops.

#include <batchlas.hh>
using namespace batchlas;
 
Queue ctx(Device::default_device());
Matrix<float> A = /* n x n, batch items */;
Matrix<float> V = /* n x k trial vectors, same batch */;
 
// Explicit backend and scalar type; MFormat is deduced from A.
Vector<float> theta = ritz_values<Backend::CUDA, float>(ctx, A, V);
float theta_j_b = theta(j, b);   // column j of batch item b

It is declared for MatrixView, plus one owning-Matrix twin so that the partially-explicit ritz_values<B, T>(ctx, A, V) compiles (include/batchlas/blas/extensions.hh:2018-2023).

Workspace form

Vector<float> theta(k, batch);   // one value per column per batch item
 
size_t bytes = ritz_values_buffer_size(ctx, A, V, theta);
UnifiedVector<std::byte> ws(bytes);
ritz_values(ctx, A, V, theta, ws);   // backend taken from the queue
ctx.wait();                          // theta is not readable before this

template <Backend B, typename T, MatrixFormat MFormat>
Event ritz_values(Queue& ctx,
                  const MatrixView<T, MFormat>& A,
                  const MatrixView<T, MatrixFormat::Dense>& V,
                  const VectorView<typename base_type<T>::type>& ritz_vals,
                  Span<std::byte> workspace);
 
template <Backend B, typename T, MatrixFormat MFormat>
size_t ritz_values_buffer_size(Queue& ctx,
                               const MatrixView<T, MFormat>& A,
                               const MatrixView<T, MatrixFormat::Dense>& V,
                               const VectorView<typename base_type<T>::type>& ritz_vals);

Both are also reachable without the backend (BATCHLAS_DISPATCH_ON_QUEUE) and take owning Matrix arguments directly. For sparse \(A\) pass a CSR matrix, e.g. auto As = A.convert_to<MatrixFormat::CSR>();.

Parameters and layout

Parameter	Meaning
A	\(n \times n\) per batch item, Dense or CSR; must be square
V	dense \(n \times k\) per batch item, same batch count as A; columns are the trial vectors
ritz_vals	real VectorView, size \(\ge k\), same batch count; entry \(j\) of item \(b\) is ritz_vals(j, b)
workspace	at least ritz_values_buffer_size(...) bytes: \(A V\) ( \(n k \cdot \text{batch}\) scalars), a batch pointer array, and the spmm workspace for CSR

Shape mismatches fail an assert, so a release build does not diagnose them.

Supported configurations

• Scalar types: float, double, std::complex<float>, std::complex<double>.
• Formats of A: Dense, CSR.
• Backends: NETLIB, CUDA, ROCM, each only when compiled in (src/ops/ritz_values/ritz_values.cc:182-192). There is no Backend::SYCL or Backend::MKL instantiation.

Tests

tests/ritz_values_tests.cc checks a diagonal matrix with identity trial vectors and a \(10 \times 10\) tridiagonal Toeplitz matrix ( \(a = 2, b = c = -1\)) at batch 2 against \(\lambda_j = a - 2\sqrt{bc}\cos\!\big(j\pi/(n+1)\big)\), over every compiled backend and scalar type. Both are dense: the CSR path has no test (open debt; spmm_tests covers spmm only).

lanczos and syevx_lobpcg.cc compute their Ritz values internally; the in-tree callers of ritz_values are the tests and the Python binding. Not implemented: an optional residual norm \(\lVert A v - \theta v \rVert\), skipping the denominator for orthonormal \(V\), and blocked evaluation for very large \(k\).