Add opnorm implementation

Project.toml

@@ -11,14 +11,18 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [weakdeps]
 AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e"
+Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
+GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb"
 LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b"
 Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
+TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e"
 
 [extensions]
 LinearOperatorsAMDGPUExt = "AMDGPU"
+LinearOperatorsOpNormExt = ["Arpack", "TSVD", "GenericLinearAlgebra"]
 LinearOperatorsCUDAExt = "CUDA"
 LinearOperatorsChainRulesCoreExt = "ChainRulesCore"
 LinearOperatorsJLArraysExt = "JLArrays"
@@ -27,14 +31,17 @@ LinearOperatorsMetalExt = "Metal"
 
 [compat]
 AMDGPU = "2.2.0"
+Arpack = "0.5.4"
 CUDA = "5.9.6"
 ChainRulesCore = "1"
 FastClosures = "0.3"
+GenericLinearAlgebra = "0.3.19"
 JLArrays = "0.3"
 LDLFactorizations = "0.10"
 LinearAlgebra = "1.10"
 Metal = "1.9.1"
 Printf = "1.10"
 SparseArrays = "1.10"
+TSVD = "0.4.4"
 TimerOutputs = "0.5"
 julia = "1.10"

ext/LinearOperatorsOpNormExt.jl

@@ -0,0 +1,138 @@
+module LinearOperatorsOpNormExt
+
+using LinearOperators
+using Arpack
+using TSVD
+using LinearAlgebra
+using GenericLinearAlgebra
+import Arpack: eigs, svds
+
+import LinearOperators: estimate_opnorm
+
+function estimate_opnorm(B::AbstractLinearOperator; kwargs...)
+  _estimate_opnorm(B, eltype(B); kwargs...)
+end
+
+# Fallback for non-standard number types (uses TSVD)
+function _estimate_opnorm(B, ::Type{T}; kwargs...) where {T}
+  _, s, _ = tsvd(B, 1; kwargs...)
+  return s[1], true
+end
+
+# Optimized path for standard FloatingPoint/Complex types
+# Note: Arpack strictly only supports Float32, Float64, ComplexF32, and ComplexF64.
+# Any other type (like BigFloat or Float16) must use the TSVD fallback above.
+function _estimate_opnorm(
+  B,
+  ::Type{T};
+  kwargs...
+) where {T <: Union{Float32, Float64, ComplexF32, ComplexF64}}
+  if ishermitian(B)
+    return opnorm_eig(B; kwargs...)
+  else
+    return opnorm_svd(B; kwargs...)
+  end
+end
+
+function opnorm_eig(B; max_attempts::Int = 3, tiny_dense_threshold = 5, kwargs...)
+  n = size(B, 1)
+
+  # Check if we can use direct dense methods (only if B supports conversion to Matrix)
+  if n ≤ tiny_dense_threshold
+    try
+      return maximum(abs, eigen(Matrix(B)).values), true
+    catch
+      # If Matrix(B) fails or isn't efficient, fall through to iterative
+    end
+  end
+
+  nev = 1
+  ncv = max(20, 2*nev + 1)
+
+  for attempt = 1:max_attempts
+    try
+      d, nconv, _, _, _ = eigs(B; nev = nev, ncv = ncv, which = :LM, ritzvec = false, check = 1, kwargs...)
+
+      if nconv == 1
+        return abs(d[1]), true
+      end
+
+    catch e
+      if e isa Arpack.ARPACKException ||
+         occursin("ARPACK", string(e)) ||
+         occursin("AUPD", string(e))
+        if ncv >= n
+          @warn "ARPACK failed and NCV cannot be increased further." exception=e
+          rethrow(e)
+        end
+      else
+        rethrow(e)
+      end
+    end
+
+    if attempt < max_attempts
+      old_ncv = ncv
+      ncv = min(2 * ncv, n)
+      if ncv > old_ncv
+        @warn "opnorm_eig: increasing NCV from $old_ncv to $ncv and retrying."
+      else
+        break
+      end
+    end
+  end
+
+  return NaN, false
+end
+
+function opnorm_svd(B; max_attempts::Int = 3, tiny_dense_threshold = 5, kwargs...)
+  m, n = size(B)
+
+  if min(m, n) ≤ tiny_dense_threshold
+    try
+      return maximum(svd(Matrix(B)).S), true
+    catch
+      # Fallback
+    end
+  end
+
+  min_dim = min(m, n)
+
+  nsv = 1
+  ncv = 10
+
+  for attempt = 1:max_attempts
+    try
+      s, nconv, _, _, _ = svds(B; nsv = nsv, ncv = ncv, ritzvec = false, check = 1, kwargs...)
+
+      if nconv >= 1
+        return maximum(s.S), true
+      end
+
+    catch e
+      if e isa Arpack.ARPACKException ||
+         occursin("ARPACK", string(e)) ||
+         occursin("AUPD", string(e))
+        if ncv >= min_dim
+          @warn "ARPACK failed and NCV cannot be increased further." exception=e
+          rethrow(e)
+        end
+      else
+        rethrow(e)
+      end
+    end
+
+    if attempt < max_attempts
+      old_ncv = ncv
+      ncv = min(2 * ncv, min_dim)
+      if ncv > old_ncv
+        @warn "opnorm_svd: increasing NCV from $old_ncv to $ncv and retrying."
+      else
+        break
+      end
+    end
+  end
+
+  return NaN, false
+end
+
+end

src/utilities.jl

@@ -1,5 +1,5 @@
 export check_ctranspose,
-  check_hermitian, check_positive_definite, normest, solve_shifted_system!, ldiv!
+  check_hermitian, check_positive_definite, normest, solve_shifted_system!, ldiv!, estimate_opnorm
 import LinearAlgebra.ldiv!
 
 """
@@ -261,8 +261,7 @@ Solves the linear system Bx = b.
 
 - `x::AbstractVector{T}`: The modified solution vector containing the solution to the linear system.
 
-### Examples:
-
+### Example:
 ```julia
 
 # Create an L-BFGS operator
@@ -276,6 +275,7 @@ b = rand(10)
 ldiv!(x, B, b)
 
 # The vector `x` now contains the solution
+```
 """
 
 function ldiv!(
@@ -287,3 +287,28 @@ function ldiv!(
   solve_shifted_system!(x, B, b, T(0.0))
   return x
 end
+
+
+"""
+    estimate_opnorm(B::AbstractLinearOperator; kwargs...)
+
+Compute the estimate of the operator 2-norm (largest singular value) of a linear operator `B`.
+
+This method dispatches to efficient algorithms depending on the properties and size of `B`.
+If `B` wraps a small dense matrix, it may use direct LAPACK routines. For larger operators,
+it uses iterative methods (ARPACK or TSVD) to estimate the norm efficiently.
+
+**Note:** This function allocates memory. It requires `Arpack.jl` and `TSVD.jl` 
+(and `GenericLinearAlgebra.jl` for generic types) to be loaded.
+
+# Arguments
+- `B::AbstractLinearOperator`: The linear operator to analyze.
+- `kwargs...`: Optional keyword arguments passed to the underlying norm estimation routines 
+  (e.g., `max_attempts`, `tiny_dense_threshold`).
+
+# Returns
+- A tuple `(norm, success)` where:
+    - `norm` is the estimated operator 2-norm of `B` (largest singular value or eigenvalue in absolute value).
+    - `success` is a boolean indicating whether the iterative method (if used) reported successful convergence.
+"""
+function estimate_opnorm end

test/Project.toml

@@ -1,5 +1,7 @@
 [deps]
 Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
+TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e"
+GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
 JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb"
@@ -15,6 +17,8 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
 Arpack = "0.5"
+TSVD = "0.4.4"
+GenericLinearAlgebra = "0.3.19"
 ChainRulesCore = "1"
 FastClosures = "0.3"
 JLArrays = "0.3"

test/runtests.jl

@@ -13,6 +13,7 @@ include("test_kron.jl")
 include("test_callable.jl")
 include("test_deprecated.jl")
 include("test_normest.jl")
+include("test_estimate_opnorm.jl")
 include("test_diag.jl")
 include("test_chainrules.jl")
 include("test_solve_shifted_system.jl")

test/test_estimate_opnorm.jl

@@ -0,0 +1,61 @@
+using TSVD
+using Arpack
+using GenericLinearAlgebra
+
+@testset "estimate_opnorm type stability and dispatch" begin
+  @testset "Type $T" for T in [Float32, Float64, ComplexF32, ComplexF64]
+    # A) Rectangular Matrix (Forces SVD path)
+    J_mat = zeros(T, 3, 2)
+    J_mat[1, 1] = 3.0
+    J_mat[2, 2] = 1.0
+    J_op = LinearOperator(J_mat)
+
+    σ, ok = estimate_opnorm(J_op)
+    @test ok
+    @test isapprox(σ, 3.0; rtol = 1e-5)
+
+    # B) Hermitian Wrapper (Forces Eig path via dispatch)
+    H_data = rand(T, 10, 10)
+    H_data = H_data + H_data'
+    H = Hermitian(H_data)
+
+    H_op = LinearOperator(H)
+
+    exact_norm = opnorm(H) # Uses standard LinearAlgebra
+    est_norm, ok = estimate_opnorm(H_op) # Uses our extension
+
+    @test ok
+    @test isapprox(est_norm, exact_norm; rtol = 1e-5)
+
+    # C) Symmetric Wrapper (Dispatch varies by Real vs Complex)
+    S_data = rand(T, 10, 10)
+    S_data = S_data + transpose(S_data)
+    S = Symmetric(S_data)
+
+    S_op = LinearOperator(S)
+
+    exact_norm_S = opnorm(Matrix(S))
+    est_norm_S, ok_S = estimate_opnorm(S_op)
+
+    @test ok_S
+    @test isapprox(est_norm_S, exact_norm_S; rtol = 1e-5)
+
+    # C.1) Specific check: Ensure Symmetric{Complex} works 
+    # (This is NOT Hermitian, so it must successfully fall back to the SVD path)
+    if T <: Complex
+      # Ensure the operator didn't falsely flag itself as Hermitian
+      @test !ishermitian(S_op)
+    end
+  end
+
+  @testset "BigFloat (Generic)" begin
+    # -------------------------
+    B_mat = Matrix{BigFloat}([2.0 0.0; 0.0 -1.0])
+    B_op = LinearOperator(B_mat)
+
+    λ_bf, ok_bf = estimate_opnorm(B_op)
+
+    @test ok_bf
+    @test isapprox(λ_bf, BigFloat(2); rtol = 1e-12)
+  end
+end