From 9cabd70f9b4116fb9ea1940fa484a4cb6bc69290 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Mon, 24 Apr 2023 14:57:52 -0400 Subject: [PATCH 1/4] use power iteration for opnorm2 --- stdlib/LinearAlgebra/src/generic.jl | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index c66f59838e8ba..c80988ee63d3a 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -670,13 +670,31 @@ function opnorm1(A::AbstractMatrix{T}) where T return convert(Tnorm, nrm) end +# Uses power iteration to compute the maximal singular value. +# falls back to svdvals if it runs into numerics issues. function opnorm2(A::AbstractMatrix{T}) where T require_one_based_indexing(A) m,n = size(A) Tnorm = typeof(float(real(zero(T)))) if m == 0 || n == 0 return zero(Tnorm) end if m == 1 || n == 1 return norm2(A) end - return svdvals(A)[1] + # to minimize the chance of x being orthogonal to the largest eigenvector + x = rand(Tnorm, m) + tmp = similar(x) + At = A' + v = one(Tnorm) + # this will converge quickly as long as the top two eigenvalues are distinct + for i in 1:n + mul!(x, At, mul!(tmp, A, x)) + newv = norm(x) + # the numerics got very wonky + !isfinite(newv) && return first(svdvals(A)) + isapprox(v, newv; rtol=10*eps(T)) && return sqrt(newv) + v = newv + x .*= inv(newv) + end + # iteration failed to converge + return first(svdvals(A)) end function opnormInf(A::AbstractMatrix{T}) where T From 0a8b9e22a4a82a9ee6a9c8f6f768479a86eb2984 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Mon, 24 Apr 2023 15:02:15 -0400 Subject: [PATCH 2/4] minor type improvement --- stdlib/LinearAlgebra/src/generic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index c80988ee63d3a..f193ab3e67fff 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -689,7 +689,7 @@ function opnorm2(A::AbstractMatrix{T}) where T newv = norm(x) # the numerics got very wonky !isfinite(newv) && return first(svdvals(A)) - isapprox(v, newv; rtol=10*eps(T)) && return sqrt(newv) + isapprox(v, newv; rtol=10*eps(Tnorm)) && return sqrt(newv) v = newv x .*= inv(newv) end From c7989025c4a37adbdc21b4b2f925a015a62cebec Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Tue, 25 Apr 2023 11:32:57 -0400 Subject: [PATCH 3/4] Nonsquare typo and use `randn` --- stdlib/LinearAlgebra/src/generic.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index f193ab3e67fff..f638dfdb66f48 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -679,8 +679,8 @@ function opnorm2(A::AbstractMatrix{T}) where T if m == 0 || n == 0 return zero(Tnorm) end if m == 1 || n == 1 return norm2(A) end # to minimize the chance of x being orthogonal to the largest eigenvector - x = rand(Tnorm, m) - tmp = similar(x) + x = randn(Tnorm, n) + tmp = zeros(Tnorm, m) At = A' v = one(Tnorm) # this will converge quickly as long as the top two eigenvalues are distinct From d885aba46f9ac59250c3d05336c0710c10abb51f Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Wed, 3 May 2023 23:56:11 -0400 Subject: [PATCH 4/4] improvements --- stdlib/LinearAlgebra/src/generic.jl | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index f638dfdb66f48..5e0ba1224cd02 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -678,20 +678,22 @@ function opnorm2(A::AbstractMatrix{T}) where T Tnorm = typeof(float(real(zero(T)))) if m == 0 || n == 0 return zero(Tnorm) end if m == 1 || n == 1 return norm2(A) end - # to minimize the chance of x being orthogonal to the largest eigenvector - x = randn(Tnorm, n) - tmp = zeros(Tnorm, m) At = A' - v = one(Tnorm) + # to minimize the chance of x being orthogonal to the largest eigenvector + x = randn(T, n) + tmp = A * x # this will converge quickly as long as the top two eigenvalues are distinct for i in 1:n - mul!(x, At, mul!(tmp, A, x)) - newv = norm(x) - # the numerics got very wonky - !isfinite(newv) && return first(svdvals(A)) - isapprox(v, newv; rtol=10*eps(Tnorm)) && return sqrt(newv) - v = newv - x .*= inv(newv) + mul!(tmp, A, x) + v1 = norm(tmp) + !isfinite(v1) && return v1 + tmp .*= inv(v1) + mul!(x, At, tmp) + v2 = norm(x) + !isfinite(v2) && return v2 + # tune this better + isapprox(v1, v2) && return v2 + x .*= inv(v2) end # iteration failed to converge return first(svdvals(A))