Skip to content

Implement the 2-norm in opnorm for sparse matrices - #859

Open
ViralBShah wants to merge 1 commit into
mainfrom
vs/opnorm2
Open

ViralBShah wants to merge 1 commit into
mainfrom
vs/opnorm2

Conversation

@ViralBShah

Copy link
Copy Markdown
Member

Fixes #119. opnorm defaults to p = 2, which threw for every sparse matrix that is not a single row or column:

julia> opnorm(sparse([1.0 2.0 0.0; 0.0 1.0 0.0; 0.0 0.0 0.3]))
ERROR: ArgumentError: 2-norm not yet implemented for sparse matrices. Try opnorm(Array(A)) or opnorm(A, p) where p=1 or Inf.

It now returns the largest singular value from Golub-Kahan-Lanczos bidiagonalization (internal helper opnorm2est), which needs only A*v and A'*u, so each step is O(nnz) and the only storage is two vectors plus the bidiagonal. The iteration stops when the residual of the leading Ritz pair, read off the SVD of the small Bidiagonal, falls below sqrt(eps()) relative to it; the error in the value is of the order of the squared residual. The result is therefore an estimate that converges from below, unlike the dense opnorm. If max(100, 2*min(m, n)) steps pass without convergence the current estimate is returned. The start vector comes from a fixed-seed Xoshiro, so repeated calls agree exactly; a fixed vector such as ones does not work, because it lies in the null space of every graph Laplacian and gives 0 for [1 -1; -1 1]. The small SVD is O(k^2), so it runs at every step up to 32 and every 8th step after that, or when the recurrence is about to break down. The return type follows the dense method. No new keywords, methods or dependencies.

Against opnorm(Array(A)) on nightly the relative error was 1e-15 to 1e-16 for real and complex Float32/Float64, Int, Bool and Rational eltypes, both rectangular orientations, low-rank matrices, repeated and clustered singular values (1 and 1 + 1e-6), all stored zeros, and entries scaled by 1e±200. A 90000×90000 2-D Laplacian takes about 1.6 s and matches the analytic value to 2e-14. test/linalg.jl passes on nightly; the old @test_throws for p = 2 is removed.

Left out: opnorm of adjoints, transposes and views of sparse matrices still reaches LinearAlgebra's generic method (elementwise for p = 1, Inf, a MethodError from svdvals! for p = 2); that is a follow-up touching dispatch. This PR wants a sentence in docs/src/index.md saying the sparse 2-norm is an iterative estimate; not backportable, since it is new behaviour.

Written by Claude Code (Claude Fable 5.1).

🤖 Generated with Claude Code

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
@codecov

codecov Bot commented Sep 19, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 92.91%. Comparing base (d0d4f2e) to head (7543037).
⚠️ Report is 2 commits behind head on main.

Additional details and impacted files
@@            Coverage Diff             @@
##             main     #859      +/-   ##
==========================================
+ Coverage   92.89%   92.91%   +0.01%     
==========================================
  Files          12       12              
  Lines        8982     9002      +20     
==========================================
+ Hits         8344     8364      +20     
  Misses        638      638              

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@ViralBShah

Copy link
Copy Markdown
Member Author

@dkarrasch Would be great for your review on this - I don't want to merge this one without review.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Operator norm (p = 2) for sparse matrices is not implemented

1 participant