# Any speed improvements to the implementation of \`BiRank\`?
**URL:** https://discourse.julialang.org/t/any-speed-improvements-to-the-implementation-of-birank/100305
**Category:** Graphs
**Tags:** graphs
**Created:** [June 13, 2023, 11:43pm UTC](https://discourse.julialang.org/t/any-speed-improvements-to-the-implementation-of-birank/100305 "2023-06-13T23:43:29Z")
**Posts on this page:** 16
**Page:** 1
### Author:  [@danielw2904](https://discourse.julialang.org/u/danielw2904)
#### Post date: [June 13, 2023, 11:43pm UTC](https://discourse.julialang.org/t/any-speed-improvements-to-the-implementation-of-birank/100305/1 "2023-06-13T23:43:29Z")
Hi All,
I implemented the paper [BiRank: Towards Ranking on Bipartite Graphs](https://ieeexplore.ieee.org/document/7572089) in Julia (R and Python implementations already exist [here](https://github.com/BrianAronson/birankr)). I was wondering if there are any obvious performance problems I missed? I noticed large differences when the `eltype` of my Sparse `W` is `UInt32` (much slower) compared to `Float64`. Also should I try to integrate it with Graphs.jl given that it is specifically for bipartite graphs? Or create a separate package?
Thanks!
```julia
module BiRank
using LinearAlgebra
using SparseArrays
using ProgressMeter
export birank
"""
birank(W; method="BiRank",
α=0.85, β=0.85, u⁰=nothing, p⁰=nothing,
max_iter=200, tol=1.0e-10)
Implements Algorithm 1 in
BiRank: Towards Ranking on Bipartite Graphs by
Xiangnan He, Ming Gao Member, Min-Yen Kan Member, and Dingxian Wang
- W: Weighted adjacency matrix (dim: |U| x |P|)
- method: one of "HITS", "CoHITS", "BGER", "BGRM", "BiRank"
- α: damping factor for 'p'
- β: damping factor for 'u'
- u⁰: query vector for 'u'. Defaults to 1/|U| for all elements.
- p⁰: query vector for 'p'. Defaults to 1/|P| for all elements.
- max_iter: maximum number of iterations
- tol: tolerance for convergence
"""
function birank(W; method="BiRank",
α=0.85, β=0.85,
u⁰=nothing, p⁰=nothing,
max_iter=200, tol=1.0e-10)
method = lowercase(method)
Wᵀ = transpose(W)
## Weighted degrees
Dₚ_vec = sum(W; dims=1) |> vec
Dᵤ_vec = sum(W; dims=2) |> vec
## Avoid division by 0
Dₚ_vec[Dₚ_vec.==zero(eltype(Dₚ_vec))] .= one(eltype(Dₚ_vec))
Dᵤ_vec[Dᵤ_vec.==zero(eltype(Dᵤ_vec))] .= one(eltype(Dᵤ_vec))
Dₚ⁻¹ = Diagonal(one(eltype(Dₚ_vec)) ./ Dₚ_vec)
Dᵤ⁻¹ = Diagonal(one(eltype(Dᵤ_vec)) ./ Dᵤ_vec)
## Table 1 in "BiRank: Towards Ranking on Bipartite Graphs"
is_hits = method == "hits"
if is_hits
S = W
Sᵀ = Wᵀ
elseif method == "cohits"
## S = W Dₚ⁻¹
## Sᵀ = Wᵀ Dᵤ⁻¹
S = W * Dₚ⁻¹
Sᵀ = Wᵀ * Dᵤ⁻¹
elseif method == "bger"
## S = Dᵤ⁻¹ W
## Sᵀ = Dₚ⁻¹ Wᵀ
S = Dᵤ⁻¹ * W
Sᵀ = Dₚ⁻¹ * Wᵀ
elseif method == "bgrm"
## S = Dᵤ⁻¹ W Dₚ⁻¹
## Sᵀ = Dₚ⁻¹ Wᵀ Dᵤ⁻¹
S = (Dᵤ⁻¹ * W) * Dₚ⁻¹
Sᵀ = transpose(S)
elseif method == "birank"
## S = sqrt(Dᵤ)⁻¹ W sqrt(Dₚ)⁻¹
## Sᵀ = sqrt(Dₚ)⁻¹ Wᵀ sqrt(Dᵤ)⁻¹
sqrtDᵤ⁻¹ = Diagonal(one(eltype(Dᵤ_vec)) ./ sqrt.(abs.(Dᵤ_vec)))
sqrtDₚ⁻¹ = Diagonal(one(eltype(Dₚ_vec)) ./ sqrt.(abs.(Dₚ_vec)))
S = (sqrtDᵤ⁻¹ * W) * sqrtDₚ⁻¹
Sᵀ = transpose(S)
else
error("""method must be one of "HITS", "CoHITS", "BGER", "BGRM", "BiRank\"""")
end
isnothing(u⁰) && (u⁰ = fill(one(eltype(Sᵀ)) / size(W, 1), size(W, 1)))
uₗ = copy(u⁰)
isnothing(p⁰) && (p⁰ = fill(one(eltype(S)) / size(W, 2), size(W, 2)))
pₗ = copy(p⁰)
normalizer_names = Dict("hits" => "HITS", "cohits" => "CoHITS", "bger" => "BGER", "bgrm" => "BGRM", "birank" => "BiRank")
progress = ProgressUnknown("Running $(normalizer_names[method])...")
let u = copy(uₗ), p = copy(pₗ)
for i in 1:max_iter
p .= α * (Sᵀ * u) + (one(α) - α) * p⁰
is_hits && (p ./= sum(p))
u .= β * (S * p) + (one(β) - β) * u⁰
is_hits && (u ./= sum(u))
εₚ = sum(abs.(p - pₗ))
εᵤ = sum(abs.(u - uₗ))
if εₚ < tol && εᵤ < tol
ProgressMeter.next!(progress)
ProgressMeter.finish!(progress)
@info "Converged after $i iterations"
break
end
ProgressMeter.next!(progress)
uₗ .= u
pₗ .= p
if i == max_iter
ProgressMeter.finish!(progress)
@warn "Not converged after $max_iter iterations. Consider increasing max_iter."
end
end
return u, p
end
end
end
```
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [June 13, 2023, 11:49pm UTC](https://discourse.julialang.org/t/any-speed-improvements-to-the-implementation-of-birank/100305/2 "2023-06-13T23:49:26Z")
One thing is that rather than passing `method` as a string, you should use a `Symbol` (or possibly dispatch). String comparisons happen letter by letter. Symbol comparisons are a single integer compare, and if you used multiple dispatch, it would all get done at compile time.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [June 13, 2023, 11:51pm UTC](https://discourse.julialang.org/t/any-speed-improvements-to-the-implementation-of-birank/100305/3 "2023-06-13T23:51:55Z")
Also, instead of `(p ./= sum(p))` you can do `p .*= inv(sum(p))` which will replace O(n) divisions with 1 division and a bunch of multiplications.
You’re also missing a number of dots which are leading to a large number of unnecessary allocations.
---
### Author:  [@danielw2904](https://discourse.julialang.org/u/danielw2904)
#### Post date: [June 14, 2023, 12:05am UTC](https://discourse.julialang.org/t/any-speed-improvements-to-the-implementation-of-birank/100305/4 "2023-06-14T00:05:58Z")
Thank you! For multiple dispatch would you use `Val{method}` or is there another way?
---
### Author:  [@danielw2904](https://discourse.julialang.org/u/danielw2904)
#### Post date: [June 14, 2023, 12:07am UTC](https://discourse.julialang.org/t/any-speed-improvements-to-the-implementation-of-birank/100305/5 "2023-06-14T00:07:51Z")
That is very cool thanks! Should I also calculate `inv.(sqrt.(abs.(Dᵤ_vec)))` instead of `one(eltype(Dᵤ_vec)) ./ sqrt.(abs.(Dᵤ_vec))`
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [June 14, 2023, 12:08am UTC](https://discourse.julialang.org/t/any-speed-improvements-to-the-implementation-of-birank/100305/6 "2023-06-14T00:08:34Z")
One commonly used interface for this is with structs. Specifically, you would define something like
```julia
abstract type BirankMethod
struct Hits end
struct Cohits end
```
(etc). and then call it as `biranck(W, BiRank())`
---
### Author:  [@danielw2904](https://discourse.julialang.org/u/danielw2904)
#### Post date: [June 14, 2023, 12:09am UTC](https://discourse.julialang.org/t/any-speed-improvements-to-the-implementation-of-birank/100305/7 "2023-06-14T00:09:04Z")
Alright, thank you!
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [June 14, 2023, 12:09am UTC](https://discourse.julialang.org/t/any-speed-improvements-to-the-implementation-of-birank/100305/8 "2023-06-14T00:09:35Z")
yes although the ones in the loop matter more.
---
### Author:  [@danielw2904](https://discourse.julialang.org/u/danielw2904)
#### Post date: [June 14, 2023, 1:39am UTC](https://discourse.julialang.org/t/any-speed-improvements-to-the-implementation-of-birank/100305/9 "2023-06-14T01:39:19Z")
I did a significant refactor based on @Oscar_Smith’s comments, thank you! Is it worth registering this as a package?
```julia
module BiRanks
using LinearAlgebra
using SparseArrays
using ProgressMeter
export birank
export BiRank, HITS, CoHITS, BGER, BGRM
abstract type BiRankMethod end
struct BiRank <: BiRankMethod end
struct HITS <: BiRankMethod end
struct CoHITS <: BiRankMethod end
struct BGER <: BiRankMethod end
struct BGRM <: BiRankMethod end
"""
birank(W, method;
α=0.85, β=0.85, u⁰=nothing, p⁰=nothing,
max_iter=200, tol=1.0e-10)
Implements Algorithm 1 in
BiRank: Towards Ranking on Bipartite Graphs by
Xiangnan He, Ming Gao Member, Min-Yen Kan Member, and Dingxian Wang
- W: Weighted adjacency matrix (dim: |U| x |P|)
- method: one of HITS(), CoHITS(), BGER(), BGRM(), BiRank(). Default: BiRank()
- α: damping factor for 'p'
- β: damping factor for 'u'
- u⁰: query vector for 'u'. Defaults to 1/|U| for all elements.
- p⁰: query vector for 'p'. Defaults to 1/|P| for all elements.
- max_iter: maximum number of iterations
- tol: tolerance for convergence
"""
birank(W) = birank(W, BiRank())
birank(W, α::Real, β::Real) = birank(W, BiRank(); α=α, β=β)
birank(W, method::BiRankMethod, α::Real, β::Real) = birank(W, method; α=α, β=β)
birank(W, method::HITS;
α=0.85, β=0.85, u⁰=nothing, p⁰=nothing,
max_iter=200, tol=1.0e-10) = birank_inner(W, W, transpose(W), α, β, u⁰, p⁰, max_iter, tol, method)
function birank(W, method::CoHITS;
α=0.85, β=0.85,
u⁰=nothing, p⁰=nothing,
max_iter=200, tol=1.0e-10)
@assert 0 <= α <= 1
@assert 0 <= β <= 1
Wᵀ, Dₚ⁻¹, Dᵤ⁻¹ = birank_inits(W, method)
## S = W Dₚ⁻¹
## Sᵀ = Wᵀ Dᵤ⁻¹
S = W * Dₚ⁻¹
Sᵀ = Wᵀ * Dᵤ⁻¹
return birank_inner(W, S, Sᵀ, α, β, u⁰, p⁰, max_iter, tol, method)
end
function birank(W, method::BGER;
α=0.85, β=0.85,
u⁰=nothing, p⁰=nothing,
max_iter=200, tol=1.0e-10)
@assert 0 <= α <= 1
@assert 0 <= β <= 1
Wᵀ, Dₚ⁻¹, Dᵤ⁻¹ = birank_inits(W, method)
## S = Dᵤ⁻¹ W
## Sᵀ = Dₚ⁻¹ Wᵀ
S = Dᵤ⁻¹ * W
Sᵀ = Dₚ⁻¹ * Wᵀ
return birank_inner(W, S, Sᵀ, α, β, u⁰, p⁰, max_iter, tol, method)
end
function birank(W, method::BGRM;
α=0.85, β=0.85,
u⁰=nothing, p⁰=nothing,
max_iter=200, tol=1.0e-10)
@assert 0 <= α <= 1
@assert 0 <= β <= 1
_, Dₚ⁻¹, Dᵤ⁻¹ = birank_inits(W, method)
## S = Dᵤ⁻¹ W Dₚ⁻¹
## Sᵀ = Dₚ⁻¹ Wᵀ Dᵤ⁻¹
S = (Dᵤ⁻¹ * W) * Dₚ⁻¹
Sᵀ = transpose(S)
return birank_inner(W, S, Sᵀ, α, β, u⁰, p⁰, max_iter, tol, method)
end
function birank(W, method::BiRank;
α=0.85, β=0.85,
u⁰=nothing, p⁰=nothing,
max_iter=200, tol=1.0e-10)
@assert 0 <= α <= 1
@assert 0 <= β <= 1
sqrtDᵤ⁻¹, sqrtDₚ⁻¹ = birank_inits(W, method)
## S = sqrt(Dᵤ)⁻¹ W sqrt(Dₚ)⁻¹
## Sᵀ = sqrt(Dₚ)⁻¹ Wᵀ sqrt(Dᵤ)⁻¹
S = (sqrtDᵤ⁻¹ * W) * sqrtDₚ⁻¹
Sᵀ = transpose(S)
return birank_inner(W, S, Sᵀ, α, β, u⁰, p⁰, max_iter, tol, method)
end
function birank_inits_common(W)
## Weighted degrees
Dₚ_vec = sum(W; dims=1) |> vec
Dᵤ_vec = sum(W; dims=2) |> vec
## Avoid division by 0
Dₚ_vec[Dₚ_vec.==zero(eltype(Dₚ_vec))] .= one(eltype(Dₚ_vec))
Dᵤ_vec[Dᵤ_vec.==zero(eltype(Dᵤ_vec))] .= one(eltype(Dᵤ_vec))
return Dₚ_vec, Dᵤ_vec
end
transpose_ifneeded(W, method::BiRankMethod) = transpose(W)
transpose_ifneeded(W, method::BGRM) = nothing
function birank_inits(W, method::BiRankMethod)
Wᵀ = transpose_ifneeded(W, method)
Dₚ_vec, Dᵤ_vec = birank_inits_common(W)
Dₚ⁻¹ = Diagonal(inv.(Dₚ_vec))
Dᵤ⁻¹ = Diagonal(inv.(Dᵤ_vec))
return Wᵀ, Dₚ⁻¹, Dᵤ⁻¹
end
function birank_inits(W, method::BiRank)
Dₚ_vec, Dᵤ_vec = birank_inits_common(W)
sqrtDᵤ⁻¹ = Diagonal(inv.(sqrt.(abs.(Dᵤ_vec))))
sqrtDₚ⁻¹ = Diagonal(inv.(sqrt.(abs.(Dₚ_vec))))
return sqrtDᵤ⁻¹, sqrtDₚ⁻¹
end
normalize_!(v, m::BiRankMethod) = nothing
normalize_!(v, m::HITS) = v .*= inv(sum(v))
function birank_inner(W, S, Sᵀ,
α, β,
u⁰, p⁰,
max_iter, tol, method::BiRankMethod)
isnothing(u⁰) && (u⁰ = fill(one(eltype(Sᵀ)) / size(W, 1), size(W, 1)))
uₗ = copy(u⁰)
isnothing(p⁰) && (p⁰ = fill(one(eltype(S)) / size(W, 2), size(W, 2)))
pₗ = copy(p⁰)
progress = ProgressUnknown("Running $method...")
let u = copy(uₗ), p = copy(pₗ)
for i in 1:max_iter
p .= α * (Sᵀ * u) + (one(α) - α) * p⁰
normalize_!(p, method)
u .= β * (S * p) + (one(β) - β) * u⁰
normalize_!(u, method)
εₚ = sum(abs.(p - pₗ))
εᵤ = sum(abs.(u - uₗ))
if εₚ < tol && εᵤ < tol
ProgressMeter.next!(progress)
ProgressMeter.finish!(progress)
@info "Converged after $i iterations"
break
end
ProgressMeter.next!(progress)
uₗ .= u
pₗ .= p
if i == max_iter
ProgressMeter.finish!(progress)
@warn "Not converged after $max_iter iterations. Consider increasing max_iter."
end
end
return u, p
end
end
end
```
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [June 14, 2023, 2:44am UTC](https://discourse.julialang.org/t/any-speed-improvements-to-the-implementation-of-birank/100305/10 "2023-06-14T02:44:24Z")
Sure! No idea what this is useful for, but it looks good to me!
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [June 14, 2023, 5:21am UTC](https://discourse.julialang.org/t/any-speed-improvements-to-the-implementation-of-birank/100305/11 "2023-06-14T05:21:16Z")
You could also contribute it to Graphs.jl?
---
### Author:  [@danielw2904](https://discourse.julialang.org/u/danielw2904)
#### Post date: [June 14, 2023, 9:18am UTC](https://discourse.julialang.org/t/any-speed-improvements-to-the-implementation-of-birank/100305/12 "2023-06-14T09:18:54Z")
I would prefer that. Are there any facilities for bipartite graphs in Graphs.jl? Right now the adjacency matrix is |U| x |P| so each have their own index rather than a global one for all vertices.
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [June 14, 2023, 10:04am UTC](https://discourse.julialang.org/t/any-speed-improvements-to-the-implementation-of-birank/100305/13 "2023-06-14T10:04:14Z")
I am not aware of anything specific for bipartite graphs, but if you find an implementation you could contribute it to my list: [The graphs ecosystem](https://discourse.julialang.org/t/the-graphs-ecosystem/99463)
Pinging @etienne_dg cause he knows better
---
### Author:  [@etienne\_dg](https://discourse.julialang.org/u/etienne_dg)
#### Post date: [June 14, 2023, 1:36pm UTC](https://discourse.julialang.org/t/any-speed-improvements-to-the-implementation-of-birank/100305/14 "2023-06-14T13:36:05Z")
We currently have some algorithms that supports only bipartite graphs, but we do not have a specific type for bipartite graphs. Since it is not that costly to check, it is currently checked at the start, and error if not bipartite (see for example [GraphsMatching.jl/src/hungarian.jl at master · JuliaGraphs/GraphsMatching.jl · GitHub](https://github.com/JuliaGraphs/GraphsMatching.jl/blob/master/src/hungarian.jl)).
Maybe we could consider adding a Trait in the Interface for bipartite graphs in a later version.
---
### Author:  [@danielw2904](https://discourse.julialang.org/u/danielw2904)
#### Post date: [June 15, 2023, 7:12am UTC](https://discourse.julialang.org/t/any-speed-improvements-to-the-implementation-of-birank/100305/15 "2023-06-15T07:12:12Z")
I’ll have a look. I’ll have to figure out what changes for a square adjacency matrix but that shouldn’t be a problem.
---
### Author:  [@danielw2904](https://discourse.julialang.org/u/danielw2904)
#### Post date: [June 16, 2023, 7:22pm UTC](https://discourse.julialang.org/t/any-speed-improvements-to-the-implementation-of-birank/100305/16 "2023-06-16T19:22:41Z")
I optimized the code quite a bit more, especially regarding allocations. The user can now set `overwrite=true` to consume the matrix `W` to save RAM. E.g. for a matrix
`776338×864116 SparseMatrixCSC{Float64, UInt32} with 46389429 stored entries` it allocates 53.14 MiB with `overwrite=true` vs. 587.325 MiB with `overwrite=false`
```julia
module BiRanks
using LinearAlgebra
using SparseArrays
using ProgressMeter
import Profile
export birank
export BiRank, HITS, CoHITS, BGER, BGRM
abstract type BiRankMethod end
struct BiRank <: BiRankMethod end
struct HITS <: BiRankMethod end
struct CoHITS <: BiRankMethod end
struct BGER <: BiRankMethod end
struct BGRM <: BiRankMethod end
"""
birank(W, method;
α=0.85, β=0.85, u⁰=nothing, p⁰=nothing,
max_iter=200, tol=1.0e-10)
Implements Algorithm 1 in
BiRank: Towards Ranking on Bipartite Graphs by
Xiangnan He, Ming Gao Member, Min-Yen Kan Member, and Dingxian Wang
- W: Weighted adjacency matrix (dim: |U| x |P|)
- method: one of HITS(), CoHITS(), BGER(), BGRM(), BiRank(). Default: BiRank()
- α: damping factor for 'p'
- β: damping factor for 'u'
- u⁰: query vector for 'u'. Defaults to 1/|U| for all elements.
- p⁰: query vector for 'p'. Defaults to 1/|P| for all elements.
- max_iter: maximum number of iterations
- tol: tolerance for convergence
"""
birank(W) = birank(W, BiRank())
birank(W, α::Real, β::Real) = birank(W, BiRank(); α=α, β=β)
birank(W, method::BiRankMethod, α::Real, β::Real) = birank(W, method; α=α, β=β)
birank(W, method::HITS;
α=0.85, β=0.85, u⁰=nothing, p⁰=nothing,
max_iter=200, tol=1.0e-10) = birank_inner(W, W, transpose(W), α, β, u⁰, p⁰, max_iter, tol, method)
function birank(W, method::CoHITS;
α=0.85, β=0.85,
u⁰=nothing, p⁰=nothing,
max_iter=200, tol=1.0e-10, overwrite=false)
Wᵀ, Dₚ⁻¹, Dᵤ⁻¹ = birank_inits(W, method)
## S = W Dₚ⁻¹
## Sᵀ = Wᵀ Dᵤ⁻¹
if !overwrite
W = copy(W)
end
S = rmul!(W, Dₚ⁻¹)
Sᵀ = rmul!(Wᵀ, Dᵤ⁻¹)
return birank_inner(W, S, Sᵀ, α, β, u⁰, p⁰, max_iter, tol, method)
end
function birank(W, method::BGER;
α=0.85, β=0.85,
u⁰=nothing, p⁰=nothing,
max_iter=200, tol=1.0e-10, overwrite=false)
Wᵀ, Dₚ⁻¹, Dᵤ⁻¹ = birank_inits(W, method)
## S = Dᵤ⁻¹ W
## Sᵀ = Dₚ⁻¹ Wᵀ
if !overwrite
W = copy(W)
end
S = lmul!(Dᵤ⁻¹, W)
Sᵀ = lmul!(Dₚ⁻¹, Wᵀ)
return birank_inner(W, S, Sᵀ, α, β, u⁰, p⁰, max_iter, tol, method)
end
function birank(W, method::BGRM;
α=0.85, β=0.85,
u⁰=nothing, p⁰=nothing,
max_iter=200, tol=1.0e-10, overwrite=false)
_, Dₚ⁻¹, Dᵤ⁻¹ = birank_inits(W, method)
## S = Dᵤ⁻¹ W Dₚ⁻¹
## Sᵀ = Dₚ⁻¹ Wᵀ Dᵤ⁻¹
if !overwrite
W = copy(W)
end
S = rmul!(lmul!(Dᵤ⁻¹, W), Dₚ⁻¹)
Sᵀ = transpose(S)
return birank_inner(W, S, Sᵀ, α, β, u⁰, p⁰, max_iter, tol, method)
end
function birank(W, method::BiRank;
α=0.85, β=0.85,
u⁰=nothing, p⁰=nothing,
max_iter=200, tol=1.0e-10, overwrite=false)
sqrtDᵤ⁻¹, sqrtDₚ⁻¹ = birank_inits(W, method)
## S = sqrt(Dᵤ)⁻¹ W sqrt(Dₚ)⁻¹
## Sᵀ = sqrt(Dₚ)⁻¹ Wᵀ sqrt(Dᵤ)⁻¹
if !overwrite
W = copy(W)
end
#S = (sqrtDᵤ⁻¹ * W) * sqrtDₚ⁻¹
S = rmul!(lmul!(sqrtDᵤ⁻¹, W), sqrtDₚ⁻¹)
Sᵀ = transpose(S)
return birank_inner(W, S, Sᵀ, α, β, u⁰, p⁰, max_iter, tol, method)
end
function birank_inits_common(W)
@assert eltype(W) <: AbstractFloat "W must be a matrix of floats"
## Weighted degrees
Dₚ_vec = sum(W; dims=1) |> vec
Dᵤ_vec = sum(W; dims=2) |> vec
## Avoid division by 0
replace_zeros!(Dₚ_vec)
replace_zeros!(Dᵤ_vec)
return Dₚ_vec, Dᵤ_vec
end
transpose_ifneeded(W, method::BiRankMethod) = transpose(W)
transpose_ifneeded(W, method::BGRM) = nothing
function birank_inits(W, method::BiRankMethod)
Wᵀ = transpose_ifneeded(W, method)
Dₚ_vec, Dᵤ_vec = birank_inits_common(W)
Dₚ⁻¹ = Diagonal(elem_inv!(Dₚ_vec))
Dᵤ⁻¹ = Diagonal(elem_inv!(Dᵤ_vec))
return Wᵀ, Dₚ⁻¹, Dᵤ⁻¹
end
function birank_inits(W, method::BiRank)
Dₚ_vec, Dᵤ_vec = birank_inits_common(W)
sqrtDᵤ⁻¹ = Diagonal(elem_sqrt_inv!(Dᵤ_vec))
sqrtDₚ⁻¹ = Diagonal(elem_sqrt_inv!(Dₚ_vec))
return sqrtDᵤ⁻¹, sqrtDₚ⁻¹
end
normalize_!(v, m::BiRankMethod) = nothing
normalize_!(v, m::HITS) = v .*= inv(sum(v))
function birank_inner(W, S, Sᵀ,
α, β,
u⁰, p⁰,
max_iter, tol, method::BiRankMethod)
@assert 0 <= α <= 1
@assert 0 <= β <= 1
isnothing(u⁰) && (u⁰ = fill(one(eltype(Sᵀ)) / size(W, 1), size(W, 1)))
uₗ = copy(u⁰)
isnothing(p⁰) && (p⁰ = fill(one(eltype(S)) / size(W, 2), size(W, 2)))
pₗ = copy(p⁰)
progress = ProgressUnknown("Running $method...")
let u = copy(uₗ), p = copy(pₗ)
for i in 1:max_iter
p .= α .* mul!(p, Sᵀ, u) .+ (one(α) - α) .* p⁰
normalize_!(p, method)
u .= β .* mul!(u, S, p) .+ (one(β) - β) .* u⁰
normalize_!(u, method)
εₚ = sum_absdiff(p, pₗ)
εᵤ = sum_absdiff(u, uₗ)
if εₚ < tol && εᵤ < tol
ProgressMeter.next!(progress)
ProgressMeter.finish!(progress)
@info "Converged after $i iterations"
break
end
ProgressMeter.next!(progress; showvalues=[("εₚ", εₚ), ("εᵤ", εᵤ)])
uₗ .= u
pₗ .= p
if i == max_iter
ProgressMeter.finish!(progress)
@warn "Not converged after $max_iter iterations. Consider increasing max_iter."
end
end
return u, p
end
end
function sum_absdiff(x::Vector{T}, y::Vector{T}) where {T}
s = zero(T)
@simd for i in eachindex(x, y)
s += abs(x[i] - y[i])
end
return s
end
function replace_zeros!(x::Vector{T}) where {T}
@simd for i in eachindex(x)
if x[i] == zero(T)
x[i] = one(T)
end
end
return x
end
function elem_inv!(x::Vector{T}) where {T}
@simd for i in eachindex(x)
x[i] = inv(x[i])
end
return x
end
function elem_sqrt_inv!(x::Vector{T}) where {T}
@simd for i in eachindex(x)
x[i] = inv(sqrt(abs(x[i])))
end
return x
end
end
```