User Guide

This guide provides detailed information about using sparse-kappa effectively.

Understanding Condition Numbers

What is a Condition Number?

The condition number of a matrix \(A\) measures how sensitive the solution of the linear system \(Ax = b\) is to perturbations in \(A\) or \(b\).

For an \(n \times n\) matrix \(A\):

  • 2-norm condition number: \(\kappa_2(A) = \|A\|_2 \cdot \|A^{-1}\|_2 = \frac{\sigma_{\max}}{\sigma_{\min}}\)

  • 1-norm condition number: \(\kappa_1(A) = \|A\|_1 \cdot \|A^{-1}\|_1\)

where \(\sigma_{\max}\) and \(\sigma_{\min}\) are the largest and smallest singular values of \(A\).

Interpretation:

  • \(\kappa(A) = 1\): Perfectly conditioned (e.g., identity matrix)

  • \(\kappa(A) < 100\): Well-conditioned

  • \(100 < \kappa(A) < 10^4\): Moderately conditioned

  • \(\kappa(A) > 10^4\): Ill-conditioned

  • \(\kappa(A) = \infty\): Singular matrix

Why Estimate Condition Numbers?

  1. Assessing numerical stability: High condition numbers indicate that numerical solutions may be unreliable

  2. Choosing solvers: Different solvers work better for different condition numbers

  3. Preconditioning: Helps design effective preconditioners

  4. Matrix analysis: Understand the properties of your linear system

Choosing Methods

Automatic Selection

The method='auto' option automatically selects the best method based on matrix properties:

cond = cond_estimate(A, method='auto')

Selection criteria:

  • Small matrices (n < 5000): Use svds (most accurate)

  • Symmetric matrices: Use eigsh

  • Very sparse matrices (density < 1%): Use golub-kahan

  • Default: Use lanczos

1-Norm Methods

Hager-Higham Algorithm

The only method for 1-norm estimation:

result = cond_estimate(A, norm=1, method='hager-higham')

Characteristics:

  • Fast iterative algorithm

  • Typically converges in 3-5 iterations

  • Based on power iteration for \(\|A^{-1}\|_1\)

  • Memory efficient

When to use:

  • When you specifically need 1-norm condition numbers

  • For sparse matrices where computing \(\|A\|_1\) is cheap

2-Norm Methods

Power Method

Simple power iteration for extreme singular values:

cond = cond_estimate(A, norm=2, method='power', max_iter=100)

Characteristics:

  • Simple and robust

  • Linear convergence

  • Good for quick rough estimates

When to use:

  • Quick estimates without high accuracy requirements

  • As a fallback when other methods fail

Lanczos Method

Uses Lanczos tridiagonalization:

cond = cond_estimate(A, norm=2, method='lanczos', num_eigenvalues=10)

Characteristics:

  • Fast convergence for well-separated eigenvalues

  • Works on \(A^T A\) (symmetric)

  • Good balance of speed and accuracy

When to use:

  • Medium-sized matrices (1k-10k)

  • When you need good accuracy

  • General purpose method

Arnoldi Method

Generalization of Lanczos for non-symmetric matrices:

cond = cond_estimate(A, norm=2, method='arnoldi', restart=20)

Characteristics:

  • More general than Lanczos

  • Higher memory requirements

  • Handles non-symmetric matrices well

When to use:

  • Non-symmetric matrices

  • When Lanczos fails to converge

Golub-Kahan Bidiagonalization

Directly computes singular values via bidiagonalization:

cond = cond_estimate(A, norm=2, method='golub-kahan', num_values=6)

Characteristics:

  • More numerically stable than forming \(A^T A\)

  • Direct SVD-based approach

  • Good for ill-conditioned matrices

When to use:

  • Large sparse matrices

  • Ill-conditioned problems

  • When numerical stability is critical

PyTorch Wrapper Methods

SVDS

Uses PyTorch’s sparse SVD:

cond = cond_estimate(A, norm=2, method='svds', num_values=10)

Characteristics:

  • Most accurate method

  • Uses optimized PyTorch implementation

  • Computes actual singular values

When to use:

  • Small to medium matrices (< 10k)

  • When you need high accuracy

  • As a reference for validating other methods

EIGSH

Uses PyTorch’s symmetric eigenvalue solver:

cond = cond_estimate(A, norm=2, method='eigsh', num_values=6)

Characteristics:

  • Optimized for symmetric matrices

  • Fast and accurate

  • Requires symmetric input

When to use:

  • Symmetric or Hermitian matrices

  • Symmetric positive definite problems

LOBPCG

Locally Optimal Block Preconditioned Conjugate Gradient:

cond = cond_estimate(A, norm=2, method='lobpcg', num_values=2)

Characteristics:

  • Efficient for a few extreme eigenvalues

  • Can use preconditioners

  • Good for very large matrices

When to use:

  • Very large matrices (> 50k)

  • When you can provide a good preconditioner

  • Computing only a few eigenvalues

Parameter Tuning

Maximum Iterations

Control convergence behavior:

# More iterations for better accuracy
cond = cond_estimate(A, max_iter=200, tol=1e-8)

# Fewer iterations for speed
cond = cond_estimate(A, max_iter=20, tol=1e-4)

Guidelines:

  • Default (100) works for most cases

  • Increase for ill-conditioned matrices

  • Decrease for rough estimates

Tolerance

Convergence tolerance:

cond = cond_estimate(A, tol=1e-6)  # Default
cond = cond_estimate(A, tol=1e-8)  # Tighter convergence

Guidelines:

  • 1e-6: Good balance (default)

  • 1e-8: High accuracy requirements

  • 1e-4: Fast approximate estimates

Number of Values

Number of singular/eigenvalues to compute:

cond = cond_estimate(A, method='svds', num_values=10)

Guidelines:

  • More values = better accuracy but slower

  • Default (6) usually sufficient

  • For large gaps between eigenvalues, fewer values may suffice

Best Practices

Performance Optimization

  1. Choose the right method for your matrix type

  2. Start with ``auto`` method selection

  3. Use symmetric methods (eigsh) when applicable

  4. Reduce ``num_values`` for faster estimates

  5. Monitor GPU memory for very large matrices

Accuracy Considerations

  1. Validate with ``svds`` on small test cases

  2. Check convergence status in results

  3. Use tight tolerances for critical applications

  4. Be aware of numerical precision limits

Memory Management

from sparse_kappa.backend import torch_api as cp

# Clear GPU memory periodically
cp.get_default_memory_pool().free_all_blocks()

# Check GPU memory usage
mempool = cp.get_default_memory_pool()
print(f"Used bytes: {mempool.used_bytes()}")
print(f"Total bytes: {mempool.total_bytes()}")

Error Handling

try:
    result = cond_estimate(A, norm=2, method='lanczos', verbose=True)
    if not result['converged']:
        print("Warning: did not converge")
except Exception as e:
    print(f"Error: {e}")
    # Fall back to different method
    result = cond_estimate(A, norm=2, method='power')

Common Pitfalls

  1. Singular matrices: Will return inf or very large numbers

  2. Memory errors: Reduce num_values or use more memory-efficient methods

  3. Non-convergence: Increase max_iter or try different method

  4. CUDA version mismatch: Ensure PyTorch matches your CUDA version

  5. Dense matrices: Convert to sparse format first for efficiency

Advanced Usage

Custom Linear Operators

from sparse_kappa.backend.sparse.linalg import LinearOperator

def matvec(v):
    # Custom matrix-vector product
    return custom_operation(v)

A_op = LinearOperator((n, n), matvec=matvec, dtype=cp.float64)
# Note: Not all methods support LinearOperators

Batch Processing

matrices = [sp.random(1000, 1000, density=0.01) for _ in range(10)]

results = []
for A in matrices:
    cond = cond_estimate(A, method='golub-kahan')
    results.append(cond)

print(f"Mean condition number: {np.mean(results):.2e}")
print(f"Max condition number: {np.max(results):.2e}")

Integration with Other Libraries

import scipy.sparse as scipy_sp

# Convert from SciPy sparse to PyTorch sparse
A_scipy = scipy_sp.random(1000, 1000, density=0.01, format='csr')
A_torch = sp.csr_matrix(cp.asarray(A_scipy.toarray()))

cond = cond_estimate(A_torch)

# Convert back if needed
result_cpu = cp.asnumpy(result['singular_values'])