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 24 stored entries:
  ⋅         ⋅        0.191601  0.41821   0.359767
 0.222009  0.467676   ⋅        0.167912  0.347202
 0.549691  0.909353  0.803885   ⋅        0.772499
  ⋅         ⋅        0.828975  0.987188   ⋅ 
  ⋅        0.828249  0.29727   0.81902   0.507742
 0.860259   ⋅        0.5879     ⋅        0.361506
 0.885582   ⋅        0.466837  0.135045  0.390669

Display of the scaled matrix A_scaled:

A_scaled

7×5 SparseArrays.SparseMatrixCSC{Float64, Int64} with 24 stored entries:
  ⋅         ⋅        0.512453  0.949043  1.0
 0.480996  0.999919   ⋅        0.344561  0.87268
 0.612546  1.0       0.999981   ⋅        0.998663
  ⋅         ⋅        0.989705  1.0        ⋅ 
  ⋅        0.999989  0.40599   0.949062  0.720661
 0.999996   ⋅        0.762871   ⋅        0.487514
 1.0        ⋅        0.588458  0.144432  0.511779

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}}:
 0.0       0.0       0.0       0.109324   0.7919     0.0        0.0
 0.0       1.01896   1.23837   0.723549   0.0        0.0        1.15952
 0.0       1.23837   1.22586   0.0        0.289815   0.0        0.0
 0.109324  0.723549  0.0       0.0        0.0        0.0        0.0775407
 0.7919    0.0       0.289815  0.0        0.0        0.0989077  0.0
 0.0       0.0       0.0       0.0        0.0989077  0.0        0.0
 0.0       1.15952   0.0       0.0775407  0.0        0.0        0.126391

Display of the scaled matrix Q_scaled:

Q_scaled

7×7 LinearAlgebra.Symmetric{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}}:
 0.0       0.0       0.0       0.188943  1.0       0.0       0.0
 0.0       0.822822  1.0       0.999984  0.0       0.0       0.999998
 0.0       1.0       0.989899  0.0       0.292658  0.0       0.0
 0.188943  0.999984  0.0       0.0       0.0       0.0       0.114452
 1.0       0.0       0.292658  0.0       0.0       0.999937  0.0
 0.0       0.0       0.0       0.0       0.999937  0.0       0.0
 0.0       0.999998  0.0       0.114452  0.0       0.0       0.116414

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