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.

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

Adding specialized ≈ and norm for structured matrices. #29724

Open
wants to merge 6 commits into
base: master
Choose a base branch
from

Conversation

mcognetta
Copy link
Contributor

@mcognetta mcognetta commented Oct 19, 2018

the == half of this PR was moved to #30108.

The logic for isapprox here is not correct yet. The way I plan to proceed is to implement specialized norm for structured matrices and then the fallback to the generic isapprox can be made much faster even if it does not fall back to element wise comparison.


Currently, == and fall back to generic methods. This PR adds special methods for structured matrix types. It also fixes an error in bidiag.jl where == fails:

julia> x = rand(3)
3-element Array{Float64,1}:
 0.7279081769922535
 0.8949580887222155
 0.7781661021218007

julia> Bu = Bidiagonal(x, zeros(2), 'U')
3×3 Bidiagonal{Float64,Array{Float64,1}}:
 0.727908  0.0        ⋅      
  ⋅        0.894958  0.0     
  ⋅         ⋅        0.778166

julia> Bl = Bidiagonal(x, zeros(2), 'L')
3×3 Bidiagonal{Float64,Array{Float64,1}}:
 0.727908   ⋅         ⋅      
 0.0       0.894958   ⋅      
  ⋅        0.0       0.778166

julia> Bu == Bl
false

julia> Matrix(Bu) == Matrix(Bl)
true

This PR

julia> x = rand(3)
3-element Array{Float64,1}:
 0.3172029052612273 
 0.8635761133663176 
 0.17013855660697486

julia> Bu = Bidiagonal(x, zeros(2), 'U')
3×3 Bidiagonal{Float64,Array{Float64,1}}:
 0.317203  0.0        ⋅      
  ⋅        0.863576  0.0     
  ⋅         ⋅        0.170139

julia> Bl = Bidiagonal(x, zeros(2), 'L')
3×3 Bidiagonal{Float64,Array{Float64,1}}:
 0.317203   ⋅         ⋅      
 0.0       0.863576   ⋅      
  ⋅        0.0       0.170139

julia> Bu == Bl
true

julia> Matrix(Bu) == Matrix(Bl)
true

Some timings (these are pessimal cases):
v1.0

julia> D=Diagonal(zeros(1000));B=Bidiagonal(zeros(1000), zeros(999), 'U');T=Tridiagonal(zeros(999), zeros(1000), zeros(999))

julia> @time D == B
  0.636254 seconds (4 allocations: 160 bytes)
true

julia> @time D == T
  0.008160 seconds (4 allocations: 160 bytes)
true

julia> @time B == T
  0.629008 seconds (4 allocations: 160 bytes)
true

This PR

julia> D=Diagonal(zeros(1000));B=Bidiagonal(zeros(1000), zeros(999), 'U');T=Tridiagonal(zeros(999), zeros(1000), zeros(999))

julia> @time D == B
  0.000007 seconds (4 allocations: 160 bytes)
true

julia> @time D == T
  0.000005 seconds (4 allocations: 160 bytes)
true

julia> @time B == T
  0.000012 seconds (4 allocations: 160 bytes)
true

(Forgot to time ):

v1.0

julia> @time D ≈ B
  1.253851 seconds (10 allocations: 24.047 KiB)
true

julia> @time D ≈ T
  0.029518 seconds (12 allocations: 39.953 KiB)
true

julia> @time B ≈ T
  0.669546 seconds (11 allocations: 32.016 KiB)
true

This PR

julia> @time D ≈ B
  0.000044 seconds (6 allocations: 8.109 KiB)
true

julia> @time D ≈ T
  0.000045 seconds (6 allocations: 8.109 KiB)
true

julia> @time B ≈ T
  0.000036 seconds (7 allocations: 16.047 KiB)
true

rtol::Real=Base.rtoldefault(promote_leaf_eltypes(A),promote_leaf_eltypes(B),atol),
nans::Bool=false, norm::Function=norm)

return iszero(B.ev) && isapprox(A.diag, B.dv; rtol=rtol, atol=atol, nans=nans)
Copy link
Contributor

Choose a reason for hiding this comment

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

Should this be a test for approximately zero when atol!=0?

Copy link
Contributor

Choose a reason for hiding this comment

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

The norm argument is unused

Copy link
Contributor

Choose a reason for hiding this comment

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

For vecnorm it would be hypot(vecnorm(B.ev), vecnorm(A.diag .- B.dv)) <= ..., but to deal with the general case I guess you need a fallback.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

vecnorm was deprecated with v0.7. The problem with using a recursive isapprox fallback is that isapprox doesn't have a kwarg for norm when the eltype is Number.

@simonbyrne simonbyrne added the linear algebra Linear algebra label Oct 31, 2018
@mcognetta mcognetta changed the title Adding specialized == and ≈ for structured matrices. Adding specialized ≈ and norm for structured matrices. Nov 21, 2018
@andreasnoack andreasnoack self-requested a review December 11, 2018 12:47
@andreasnoack
Copy link
Member

This needs a rebase.

@KristofferC KristofferC mentioned this pull request Dec 12, 2018
52 tasks
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants