From eaac7fcba1096428566e0228fa7038d80f92a549 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 18 Jan 2022 14:28:12 +0100 Subject: [PATCH 1/9] add _orth and rewrite image_basis! for fp matrices --- src/matrix_projections.jl | 49 ++++++++++++++++++++++++++++++++------- 1 file changed, 40 insertions(+), 9 deletions(-) diff --git a/src/matrix_projections.jl b/src/matrix_projections.jl index d8dd0fb..5a51810 100644 --- a/src/matrix_projections.jl +++ b/src/matrix_projections.jl @@ -110,7 +110,6 @@ end _mproj_fitsT!(args...) = _mproj_outsT!(args...) - function _mproj_outsT!(mproj::AbstractMatrix{T}, χ::Character) where {T} mproj .= sum( c .* matrix_projection_irr(Character(table(χ), i)) for @@ -337,8 +336,11 @@ For characters or algebra elements return a basis of the row-space of `matrix_projection([hom, ]χ)` or `matrix_representation([hom, ]α)`. Two methods are employed to achieve the goal. -* By default (symbolic, exact) row echelon form is produced, and therefore there is no guarantee on the orthogonality of the returned basis vectors (rows). -* If `eltype(A) <: AbstractFloat` this is followed by a call to `svd` and the appropriate rows of its `Vt` are returned, thus the basis is (numerically) orthonormal. +* By default (symbolic, exact) row echelon form is produced, and therefore there +is no guarantee on the orthogonality of the returned basis vectors (rows). +* If `eltype(A) <: AbstractFloat` this is followed by a call to `svd` (or +`qr` in the sparse case) and the appropriate rows of its orthogonal factor are +returned, thus the basis is (numerically) orthonormal. # Examples: ```julia @@ -396,12 +398,41 @@ function image_basis!(A::AbstractSparseMatrix) return A, p end -function image_basis!(A::AbstractMatrix{T}) where {T<:Union{AbstractFloat,Complex}} +function _orth!(M::AbstractMatrix{T}) where T<:Union{AbstractFloat, Complex} + F = svd!(convert(Matrix{T}, M)) + M_rank = let m = T <: Complex ? 2eps(real(T)) : eps(T) + sum(F.S .> maximum(size(M)) * m) + end + return F.Vt, M_rank +end + +function _orth!(M::AbstractSparseMatrix{T}) where T<:Union{AbstractFloat, Complex} + F = qr(M) + M_rank= rank(F) + result = let tmp = F.Q * Matrix(I, size(F.Q, 2), M_rank) + sparse(tmp[invperm(F.prow), :]') + end + result = let m = T <: Complex ? 2eps(real(T)) : eps(T) + droptol!(result, maximum(size(result)) * m) + end + return result, M_rank +end + +function image_basis!( + A::AbstractMatrix{T}, +) where {T<:Union{AbstractFloat,Complex}} A, p = row_echelon_form!(A) - fact = svd!(convert(Matrix{T}, @view(A[1:length(p), :]))) - m = T <: Complex ? 2eps(real(T)) : eps(T) - A_rank = sum(fact.S .> maximum(size(A)) * m) - @assert A_rank == length(p) - return fact.Vt, 1:length(p) + A_orth, A_rank = _orth!(@view A[1:length(p), :]) + @assert A_rank == length(p) "_orth rank doesn't agree with rref rank!" + return A_orth, 1:A_rank end +function image_basis!( + A::AbstractSparseMatrix{T}, +) where {T<:Union{AbstractFloat,Complex}} + droptol!(A, eps(real(eltype(A))) * size(A, 1)) + # A, p = row_echelon_form!(A) + A_orth, A_rank = _orth!(A) + # @assert A_rank == length(p) "_orth rank doesn't agree with rref rank!" + return A_orth, 1:A_rank +end From 6d294aabc27ccecb7c52004cdab3ba63b9781f7c Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 18 Jan 2022 14:29:45 +0100 Subject: [PATCH 2/9] make starting_at a kwarg --- src/Characters/echelon_form.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/Characters/echelon_form.jl b/src/Characters/echelon_form.jl index 1aa4b6f..b4276c6 100644 --- a/src/Characters/echelon_form.jl +++ b/src/Characters/echelon_form.jl @@ -13,7 +13,7 @@ end Base.@propagate_inbounds function _mul_row!( A::AbstractMatrix, val, - row_idx, + row_idx; starting_at = 1, ) @boundscheck @assert 1 ≤ starting_at ≤ size(A, 2) @@ -28,7 +28,7 @@ end Base.@propagate_inbounds function _find_pivot( A::AbstractMatrix, - col_idx, + col_idx; starting_at = 1, ) k = findnext(!iszero, @view(A[:, col_idx]), starting_at) @@ -39,7 +39,7 @@ end Base.@propagate_inbounds function _reduce_column_by_pivot!( A::AbstractMatrix, row_idx, - col_idx, + col_idx; starting_at = 1, ) @boundscheck checkbounds(A, row_idx, col_idx) @@ -84,8 +84,8 @@ end Base.@propagate_inbounds function _find_pivot( A::AbstractMatrix{T}, - col_idx, - starting_at = 1; + col_idx; + starting_at = 1, atol = eps(real(eltype(A))) * size(A, 1), ) where {T<:Union{AbstractFloat,Complex{<:AbstractFloat}}} isempty(starting_at:size(A, 1)) && return false, starting_at @@ -103,7 +103,7 @@ end Base.@propagate_inbounds function _reduce_column_by_pivot!( A::AbstractMatrix{T}, row_idx, - col_idx, + col_idx; starting_at = 1, atol = eps(real(eltype(A))) * size(A, 1), ) where {T<:Union{AbstractFloat,Complex{<:AbstractFloat}}} @@ -129,7 +129,7 @@ Base.@propagate_inbounds function row_echelon_form!(A::AbstractMatrix) pivots = Int[] row_idx = 0 # also: current_rank @inbounds for col_idx in 1:size(A, 2) - found, j = _find_pivot(A, col_idx, row_idx + 1) + found, j = _find_pivot(A, col_idx, starting_at = row_idx + 1) found || continue row_idx += 1 push!(pivots, col_idx) @@ -141,12 +141,12 @@ Base.@propagate_inbounds function row_echelon_form!(A::AbstractMatrix) w = inv(A[row_idx, col_idx]) # multiply A[row_idx, :] by w - _mul_row!(A, w, row_idx, col_idx) + _mul_row!(A, w, row_idx, starting_at = col_idx) # A[current_rank, col_idx] is 1 now # zero the whole col_idx-th column above and below pivot: # to the left of col_idx everything is zero - _reduce_column_by_pivot!(A, row_idx, col_idx, col_idx) + _reduce_column_by_pivot!(A, row_idx, col_idx, starting_at = col_idx) end return A, pivots end From 82944beca02c2d1e2cb4544bb0277d7027b7489a Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 18 Jan 2022 14:30:43 +0100 Subject: [PATCH 3/9] add const FloarOrComplex union --- src/Characters/echelon_form.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/Characters/echelon_form.jl b/src/Characters/echelon_form.jl index b4276c6..ffb2f25 100644 --- a/src/Characters/echelon_form.jl +++ b/src/Characters/echelon_form.jl @@ -81,13 +81,14 @@ Base.@propagate_inbounds function __findmax( return midx, mval end +const FloatOrComplex = Union{AbstractFloat,Complex{<:AbstractFloat}} Base.@propagate_inbounds function _find_pivot( A::AbstractMatrix{T}, col_idx; starting_at = 1, atol = eps(real(eltype(A))) * size(A, 1), -) where {T<:Union{AbstractFloat,Complex{<:AbstractFloat}}} +) where {T<:FloatOrComplex} isempty(starting_at:size(A, 1)) && return false, starting_at @boundscheck @assert 1 ≤ starting_at # find the largest entry in the column below the last pivot @@ -106,7 +107,7 @@ Base.@propagate_inbounds function _reduce_column_by_pivot!( col_idx; starting_at = 1, atol = eps(real(eltype(A))) * size(A, 1), -) where {T<:Union{AbstractFloat,Complex{<:AbstractFloat}}} +) where {T<:FloatOrComplex} @boundscheck checkbounds(A, row_idx, col_idx) @boundscheck 1 ≤ starting_at ≤ size(A, 2) From 60e0b1f700a720920e679c1e71ecb74c9b5c5f58 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 18 Jan 2022 14:31:31 +0100 Subject: [PATCH 4/9] remove own __findmax --- src/Characters/echelon_form.jl | 26 ++------------------------ 1 file changed, 2 insertions(+), 24 deletions(-) diff --git a/src/Characters/echelon_form.jl b/src/Characters/echelon_form.jl index ffb2f25..a5fb93d 100644 --- a/src/Characters/echelon_form.jl +++ b/src/Characters/echelon_form.jl @@ -56,31 +56,9 @@ Base.@propagate_inbounds function _reduce_column_by_pivot!( return A end -# version over AbstractFloat -Base.@propagate_inbounds function __findmax( - f, - A::AbstractArray{T}; - atol=eps(real(eltype(A)))*max(size(A)...) -) where {T<:Union{AbstractFloat,Complex{<:AbstractFloat}}} - @assert !isempty(A) - - midx = first(eachindex(A)) - mval = f(A[midx]) - - @inbounds for idx in eachindex(A) - iszero(A[idx]) && continue - v = f(A[idx]) - if v > mval - midx, mval = idx, v - end - if v < atol - A[idx] = zero(T) - end - end +# version over AbstractFloat - return midx, mval -end const FloatOrComplex = Union{AbstractFloat,Complex{<:AbstractFloat}} Base.@propagate_inbounds function _find_pivot( @@ -94,7 +72,7 @@ Base.@propagate_inbounds function _find_pivot( # find the largest entry in the column below the last pivot @boundscheck @assert 1 ≤ col_idx ≤ size(A, 2) - midx, mval = __findmax(abs, @view(A[starting_at:end, col_idx]), atol=atol) + mval, midx = findmax(abs, @view A[starting_at:end, col_idx]) if mval < atol # the largest entry is below threshold so we zero everything in the column return false, starting_at end From c8c6fae4d311971c240cabb44ee5cd8def557009 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 18 Jan 2022 14:33:32 +0100 Subject: [PATCH 5/9] add additional _finalize_pivot_reduce! it seems that droptol! after each reduce is much faster for sparse matrices --- src/Characters/echelon_form.jl | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/Characters/echelon_form.jl b/src/Characters/echelon_form.jl index a5fb93d..acba9dd 100644 --- a/src/Characters/echelon_form.jl +++ b/src/Characters/echelon_form.jl @@ -56,6 +56,7 @@ Base.@propagate_inbounds function _reduce_column_by_pivot!( return A end +_finalize_pivot_reduce!(A::AbstractMatrix, pivot) = A # version over AbstractFloat @@ -79,6 +80,15 @@ Base.@propagate_inbounds function _find_pivot( return true, oftype(starting_at, starting_at + midx - 1) end +function _finalize_pivot_reduce!( + A::AbstractSparseMatrix{T}, + pivot, +) where {T<:FloatOrComplex} + m = T <: Complex ? 2abs(eps(T)) : eps(T) + droptol!(A, max(size(A)...) * m) + return A +end + Base.@propagate_inbounds function _reduce_column_by_pivot!( A::AbstractMatrix{T}, row_idx, @@ -126,6 +136,7 @@ Base.@propagate_inbounds function row_echelon_form!(A::AbstractMatrix) # zero the whole col_idx-th column above and below pivot: # to the left of col_idx everything is zero _reduce_column_by_pivot!(A, row_idx, col_idx, starting_at = col_idx) + A = _finalize_pivot_reduce!(A, col_idx) end return A, pivots end From d6ec89a6f3c9ff878f674b93c9c801ec2be1e054 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 18 Jan 2022 14:34:04 +0100 Subject: [PATCH 6/9] cosmetics --- src/Characters/echelon_form.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/Characters/echelon_form.jl b/src/Characters/echelon_form.jl index acba9dd..d0ed264 100644 --- a/src/Characters/echelon_form.jl +++ b/src/Characters/echelon_form.jl @@ -75,6 +75,7 @@ Base.@propagate_inbounds function _find_pivot( mval, midx = findmax(abs, @view A[starting_at:end, col_idx]) if mval < atol # the largest entry is below threshold so we zero everything in the column + @view(A[starting_at:end, col_idx]) .= zero(T) return false, starting_at end return true, oftype(starting_at, starting_at + midx - 1) @@ -106,9 +107,8 @@ Base.@propagate_inbounds function _reduce_column_by_pivot!( @inbounds A[ridx, col_idx] = zero(T) continue end - for cidx in starting_at:size(A, 2) - @inbounds A[ridx, cidx] -= - cidx == col_idx ? A[ridx, cidx] : v * A[row_idx, cidx] + @inbounds for cidx in starting_at:size(A, 2) + A[ridx, cidx] -= v * A[row_idx, cidx] end end return A @@ -121,21 +121,22 @@ Base.@propagate_inbounds function row_echelon_form!(A::AbstractMatrix) found, j = _find_pivot(A, col_idx, starting_at = row_idx + 1) found || continue row_idx += 1 + # @info col_idx push!(pivots, col_idx) # swap the rows so that A[row_idx, :] is the row with leading nz in # i-th column - _swap_rows!(A, row_idx, j) + A = _swap_rows!(A, row_idx, j) w = inv(A[row_idx, col_idx]) # multiply A[row_idx, :] by w - _mul_row!(A, w, row_idx, starting_at = col_idx) + A = _mul_row!(A, w, row_idx, starting_at = col_idx) # A[current_rank, col_idx] is 1 now # zero the whole col_idx-th column above and below pivot: # to the left of col_idx everything is zero - _reduce_column_by_pivot!(A, row_idx, col_idx, starting_at = col_idx) + A = _reduce_column_by_pivot!(A, row_idx, col_idx, starting_at = col_idx) A = _finalize_pivot_reduce!(A, col_idx) end return A, pivots From 37f2c7d450bd5d79cc092ae62bb8c9c9e8fec05a Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 18 Jan 2022 14:34:32 +0100 Subject: [PATCH 7/9] threadin in CachedExtensionHomomorphism --- src/actions.jl | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/actions.jl b/src/actions.jl index 2a9cb67..7bfdc70 100644 --- a/src/actions.jl +++ b/src/actions.jl @@ -137,8 +137,15 @@ function CachedExtensionHomomorphism( if precompute # to make sure that the future access to chom is read-only, i.e. thread-safe # one may choose to precompute the images - for g in G - induce(chom, g) + + lck = Threads.SpinLock() + @sync for g in G + Threads.@spawn begin + val = induce(action, chom.ehom, g) + lock(lck) do + chom.cache[g] = val + end + end end end return chom From 39eeb40dce9c9324b1b96dd9a23bc3e329b49cb6 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 18 Jan 2022 15:13:21 +0100 Subject: [PATCH 8/9] relax scs less on robinson form --- test/action_dihedral.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/action_dihedral.jl b/test/action_dihedral.jl index 9314148..dab92dc 100644 --- a/test/action_dihedral.jl +++ b/test/action_dihedral.jl @@ -15,7 +15,7 @@ const OPTIMIZER = optimizer_with_attributes( SCS.Optimizer, "acceleration_lookback" => 10, "max_iters" => 3_000, - "alpha" => 1.2, + "alpha" => 1.0, "eps" => 1e-6, "linear_solver" => SCS.DirectSolver, ) From 3ff20082919abc209985cb43a7d1f74ff0050c32 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 18 Jan 2022 15:18:51 +0100 Subject: [PATCH 9/9] add _eps --- src/matrix_projections.jl | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/src/matrix_projections.jl b/src/matrix_projections.jl index 5a51810..7816f3b 100644 --- a/src/matrix_projections.jl +++ b/src/matrix_projections.jl @@ -398,23 +398,23 @@ function image_basis!(A::AbstractSparseMatrix) return A, p end -function _orth!(M::AbstractMatrix{T}) where T<:Union{AbstractFloat, Complex} +_eps(::Type{T}) where {T} = T <: Complex ? 2eps(real(T)) : eps(T) + +function _orth!(M::AbstractMatrix{T}) where {T<:Union{AbstractFloat,Complex}} F = svd!(convert(Matrix{T}, M)) - M_rank = let m = T <: Complex ? 2eps(real(T)) : eps(T) - sum(F.S .> maximum(size(M)) * m) - end + M_rank = count(>(maximum(size(M)) * _eps(T)), F.S) return F.Vt, M_rank end -function _orth!(M::AbstractSparseMatrix{T}) where T<:Union{AbstractFloat, Complex} +function _orth!( + M::AbstractSparseMatrix{T}, +) where {T<:Union{AbstractFloat,Complex}} F = qr(M) - M_rank= rank(F) + M_rank = rank(F) result = let tmp = F.Q * Matrix(I, size(F.Q, 2), M_rank) sparse(tmp[invperm(F.prow), :]') end - result = let m = T <: Complex ? 2eps(real(T)) : eps(T) - droptol!(result, maximum(size(result)) * m) - end + result = droptol!(result, maximum(size(result)) * _eps(T)) return result, M_rank end @@ -430,7 +430,8 @@ end function image_basis!( A::AbstractSparseMatrix{T}, ) where {T<:Union{AbstractFloat,Complex}} - droptol!(A, eps(real(eltype(A))) * size(A, 1)) + N = LinearAlgebra.checksquare(A) + droptol!(A, N * _eps(T)) # A, p = row_echelon_form!(A) A_orth, A_rank = _orth!(A) # @assert A_rank == length(p) "_orth rank doesn't agree with rref rank!"