Equilibration algorithm

The equilibration algorithm scales a matrix so that its rows and columns have an infinite norm of 1.

using LinearAlgebra, SparseArrays
using OperatorScaling
T = Float64
ϵ = 1.0e-4 # tolerance
m, n = 7, 5
A = sprand(T, m, n, 0.6)
A_scaled, D1, D2 = equilibrate(A, ϵ = ϵ)
norm(D1 * A * D2 - A_scaled) ≤ sqrt(eps(T)) * norm(A)

true

Display of the input matrix A:

A

7×5 SparseArrays.SparseMatrixCSC{Float64, Int64} with 27 stored entries:
 0.138513  0.914493   0.15369   0.997953  0.945509
 0.900291  0.0859649   ⋅        0.791852  0.588885
 0.562162  0.0615039   ⋅        0.175462  0.45583
 0.881811  0.83381    0.552083  0.404347   ⋅ 
 0.741393   ⋅         0.276675  0.57087   0.890707
 0.637185   ⋅         0.911378  0.604656  0.5081
 0.573784   ⋅          ⋅         ⋅        0.230005

Display of the scaled matrix A_scaled:

A_scaled

7×5 SparseArrays.SparseMatrixCSC{Float64, Int64} with 27 stored entries:
 0.146131  0.999989   0.161154  1.0       0.987301
 1.0       0.0989692   ⋅        0.835405  0.647408
 0.999943  0.113391    ⋅        0.296438  0.802504
 0.999997  0.980059   0.622257  0.435525   ⋅ 
 0.840975   ⋅         0.311923  0.615047  1.0
 0.703436   ⋅         1.0       0.634022  0.555186
 0.999945   ⋅          ⋅         ⋅        0.39673

The in-place version uses storage diagonal matrices, and updates A, D1 and D2 in-place.

D1, R_k = Diagonal(Vector{T}(undef, m)), Diagonal(Vector{T}(undef, m))
D2, C_k = Diagonal(Vector{T}(undef, n)), Diagonal(Vector{T}(undef, n))
A_scaled2 = copy(A)
equilibrate!(A_scaled2, D1, D2, R_k, C_k, ϵ = ϵ)
# A_scaled2, D1 and D2 are now updated
norm(D1 * A * D2 - A_scaled2) ≤ sqrt(eps(T)) * norm(A)

true

This packages also features an implementation for symmetric matrices:

A = sprand(Float64, m, m, 0.3)
Q = Symmetric(tril(A + A'), :L)
Q_scaled, D = equilibrate(Q, ϵ = ϵ)
norm(D * Q * D - Q_scaled) ≤ sqrt(eps(T)) * norm(Q)

true

Display of the input matrix Q:

Q

7×7 LinearAlgebra.Symmetric{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}}:
 1.15324   0.0        0.0       0.0       0.0        0.438182  0.0
 0.0       0.0        1.33698   0.0       0.0503179  0.838449  0.646924
 0.0       1.33698    0.0       0.664477  0.0        0.0       0.0
 0.0       0.0        0.664477  1.23464   0.143839   0.0       0.0
 0.0       0.0503179  0.0       0.143839  0.0        0.128318  0.669058
 0.438182  0.838449   0.0       0.0       0.128318   0.0       0.0
 0.0       0.646924   0.0       0.0       0.669058   0.0       0.550066

Display of the scaled matrix Q_scaled:

Q_scaled

7×7 LinearAlgebra.Symmetric{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}}:
 1.0       0.0        0.0       0.0       0.0        0.562672  0.0
 0.0       0.0        1.0       0.0       0.0532021  0.999943  0.684005
 0.0       1.0        0.0       0.517187  0.0        0.0       0.0
 0.0       0.0        0.517187  1.0       0.158261   0.0       0.0
 0.0       0.0532021  0.0       0.158261  0.0        0.21633   1.0
 0.562672  0.999943   0.0       0.0       0.21633    0.0       0.0
 0.0       0.684005   0.0       0.0       1.0        0.0       0.82215

Q_scaled2 = copy(Q)
# D diagonal matrix and storage diagonal matrix (same size as Q)
D, C_k = Diagonal(Vector{T}(undef, m)), Diagonal(Vector{T}(undef, m))
equilibrate!(Q_scaled2, D, C_k, ϵ = ϵ)
norm(D * Q * D - Q_scaled2) ≤ sqrt(eps(T)) * norm(Q)

true