Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

use power iteration for opnorm2 #49487

Closed
wants to merge 4 commits into from
Closed

use power iteration for opnorm2 #49487

Changes from 3 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
9cabd70
use power iteration for opnorm2
oscardssmith Apr 24, 2023
0a8b9e2
minor type improvement
oscardssmith Apr 24, 2023
c798902
Nonsquare typo and use `randn`
oscardssmith Apr 25, 2023
d885aba
improvements
oscardssmith May 4, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 19 additions & 1 deletion stdlib/LinearAlgebra/src/generic.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 = randn(Tnorm, n)
oscardssmith marked this conversation as resolved.
Show resolved Hide resolved
tmp = zeros(Tnorm, m)
At = A'
v = one(Tnorm)
# this will converge quickly as long as the top two eigenvalues are distinct
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Couldn't we use a Lanczos iteration to be more robust and converge quickly even if a few eigenvalues are similar in magnitude?

Copy link
Member

@stevengj stevengj Apr 29, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This paper seems to be a popular algorithm: https://doi.org/10.1137/04060593X and is even fancier (restarting, bidiagonalization).

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

good point. There's a python implementation of that paper with a good license here: https://github.com/bwlewis/irlbpy/blob/master/irlb/irlb.py. I'll see how it does.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

good news and bad news. Good news is I have it working. Bad news is it takes about 84 lines of code. Performance testing to come.

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))
oscardssmith marked this conversation as resolved.
Show resolved Hide resolved
isapprox(v, newv; rtol=10*eps(Tnorm)) && 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
Expand Down