Empirically verify the perturbation sensitivity of doubly stochastic matrices

Authors

Affiliation

2. Numerical Stability Considerations

Numerical stability refers to how errors—stemming from finite precision arithmetic, rounding, and truncation—propagate through computations involving the matrix.

a. Conditioning

    •   Condition Number: The condition number of a matrix measures how much the output value can change for a small change in the input. For doubly stochastic matrices, the condition number can vary widely depending on the specific matrix. Some are well-conditioned (small condition numbers), implying that they are relatively stable under numerical perturbations, while others are ill-conditioned (large condition numbers), making them sensitive to errors.
    •   Perturbation Sensitivity: Doubly stochastic matrices can be sensitive to perturbations, especially those near the boundary of the Birkhoff polytope (i.e., those close to permutation matrices). Small numerical errors might push such matrices outside the set of doubly stochastic matrices or significantly alter their properties.

4. Theoretical Insights and Research

Research into the numerical stability of doubly stochastic matrices highlights several key points:

    •   Birkhoff Polytope Structure: The geometric properties of the Birkhoff polytope influence stability. Matrices that lie deep within the polytope (far from permutation matrices) tend to be better conditioned than those near the vertices.
    •   Spectral Properties: The eigenvalues and singular values of doubly stochastic matrices play a role in their numerical behavior. For instance, the presence of eigenvalues close to zero can indicate potential instability.
    •   Random Doubly Stochastic Matrices: Probabilistic analyses suggest that randomly generated doubly stochastic matrices often exhibit better conditioning on average, though specific instances may still be problematic.

using LinearAlgebra
using Random

2. Generating a Doubly Stochastic Matrix

There are multiple ways to generate a doubly stochastic matrix. One common method is using the Sinkhorn-Knopp algorithm, which iteratively normalizes the rows and columns of a non-negative matrix to make it doubly stochastic.

function sinkhorn_knopp(A; max_iters=1000, tol=1e-9, normalize_rows=true, normalize_cols=true)
    """
    Applies the Sinkhorn-Knopp algorithm to matrix A to make it doubly stochastic.

    # Arguments
    - `A::Matrix{Float64}`: Non-negative input matrix
    - `max_iters::Int`: Maximum number of iterations
    - `tol::Float64`: Tolerance for convergence
    - `normalize_rows::Bool`: Whether to normalize rows
    - `normalize_cols::Bool`: Whether to normalize columns

    # Returns
    - `Matrix{Float64}`: Doubly stochastic matrix
    """
    A = copy(A)
    n, m = size(A)
    if n != m
        error("Matrix must be square.")
    end
    for iter in 1:max_iters
        if normalize_rows
            row_sums = sum(A, dims=2)
            A = A ./ row_sums
        end
        if normalize_cols
            col_sums = sum(A, dims=1)
            A = A ./ col_sums
        end
        # Check convergence
        row_diff = normalize_rows ? maximum(abs.(sum(A, dims=2) .- 1)) : 0.0
        col_diff = normalize_cols ? maximum(abs.(sum(A, dims=1) .- 1)) : 0.0
        if row_diff < tol && col_diff < tol
            println("Converged in $iter iterations.")
            break
        end
        if iter == max_iters
            println("Reached maximum iterations without full convergence.")
        end
    end
    return A
end

# Example: Generate a random non-negative matrix and make it doubly stochastic
Random.seed!(123)  # For reproducibility
n = 5  # Size of the matrix
A = rand(n, n)
A_ds = sinkhorn_knopp(A)
println("Doubly Stochastic Matrix A_ds:\n")
display(A_ds)
println("Row sums: ", sum(A_ds, dims=2))
println("Column sums: ", sum(A_ds, dims=1))

i = 5
A_ds_i = A_ds^i
println("Matrix A_ds^$i:\n")
display(A_ds_i)
println("Row sums: ", sum(A_ds_i, dims=2))
println("Column sums: ", sum(A_ds_i, dims=1))

3. Introducing Perturbations

To assess perturbation sensitivity, we can introduce small random perturbations to the doubly stochastic matrix and observe how its properties change.

function perturb_matrix(A::Matrix{Float64}, noise_level::Float64)
    """
    Adds Gaussian noise to matrix A and projects it back to doubly stochastic.
    
    # Arguments
    - `A::Matrix{Float64}`: Original doubly stochastic matrix
    - `noise_level::Float64`: Standard deviation of Gaussian noise
    
    # Returns
    - `Matrix{Float64}`: Perturbed doubly stochastic matrix
    """
    n, m = size(A)
    noise = randn(n, m) * noise_level
    A_perturbed = A + noise
    # Ensure non-negativity
    A_perturbed[A_perturbed .< 0] .= 0.0
    # Re-apply Sinkhorn-Knopp to make it doubly stochastic
    A_perturbed_ds = sinkhorn_knopp(A_perturbed)
    return A_perturbed_ds
end

function perturb_right_stochastic(A::Matrix{Float64}, noise_level::Float64)
    """
    Perturbs a right stochastic matrix while preserving row sums.

    # Arguments
    - `A::Matrix{Float64}`: Original right stochastic matrix
    - `noise_level::Float64`: Standard deviation of Gaussian noise

    # Returns
    - `Matrix{Float64}`: Perturbed right stochastic matrix
    """
    n, m = size(A)
    # Add Gaussian noise
    noise = randn(n, m) * noise_level
    A_perturbed = A + noise
    # Set negative entries to a small positive value to maintain non-negativity
    A_perturbed[A_perturbed .< 1e-8] .= 1e-8
    # Renormalize rows to sum to 1
    row_sums = sum(A_perturbed, dims=2)
    A_perturbed_rs = A_perturbed ./ row_sums
    return A_perturbed_rs
end


# Example usage
noise_level = 0.05
A_rs_perturbed = perturb_right_stochastic(A_ds, noise_level)
println("Perturbed Right Stochastic Matrix A_rs_perturbed:\n", A_rs_perturbed)
println("Row sums: ", sum(A_rs_perturbed, dims=2))
println("Column sums: ", sum(A_rs_perturbed, dims=1))

A_rs_perturbed_i = A_rs_perturbed^100
println("Matrix A_rs_perturbed^$i:\n")
display(A_rs_perturbed_i)
println("Row sums: ", sum(A_rs_perturbed_i, dims=2))
println("Column sums: ", sum(A_rs_perturbed_i, dims=1)) 

function compute_metrics(A_original::Matrix{Float64}, A_perturbed::Matrix{Float64})
    """
    Computes various metrics to assess the difference between two matrices.
    
    # Arguments
    - `A_original::Matrix{Float64}`: Original doubly stochastic matrix
    - `A_perturbed::Matrix{Float64}`: Perturbed doubly stochastic matrix
    
    # Returns
    - `Dict`: Dictionary containing computed metrics
    """
    difference = A_original - A_perturbed
    absolute_diff = maximum(abs.(difference))
    relative_diff = maximum(abs.(difference) ./ A_original)
    
    # Condition numbers
    cond_original = cond(A_original)
    cond_perturbed = cond(A_perturbed)
    
    # Eigenvalues
    eigen_original = eigen(A_original).values
    eigen_perturbed = eigen(A_perturbed).values
    eigen_diff = maximum(abs.(eigen_original - eigen_perturbed))
    
    # Singular values
    singular_original = svdvals(A_original)
    singular_perturbed = svdvals(A_perturbed)
    singular_diff = maximum(abs.(singular_original - singular_perturbed))
    
    # Frobenius norm
    # frob_norm = norm(difference, Frobenius)
    
    return Dict(
        :absolute_difference => absolute_diff,
        :relative_difference => relative_diff,
        :condition_number_original => cond_original,
        :condition_number_perturbed => cond_perturbed,
        :max_eigenvalue_difference => eigen_diff,
        :max_singular_difference => singular_diff
        # :frobenius_norm => frob_norm
    )
end

# Compute metrics
metrics = compute_metrics(A_ds, A_perturbed)
println("Perturbation Metrics:")
for (key, value) in metrics
    println("$key: $value")
end

using CairoMakie
# Define a range of noise levels
noise_levels = [0.001, 0.005, 0.01, 0.05, 0.1]

# Initialize containers for metrics
absolute_diffs = Float64[]
relative_diffs = Float64[]
cond_originals = Float64[]
cond_perturbed_vals = Float64[]
max_eigen_diffs = Float64[]
max_singular_diffs = Float64[]
frob_norms = Float64[]

for noise in noise_levels
    A_p = perturb_matrix(A_ds, noise)
    m = compute_metrics(A_ds, A_p)
    push!(absolute_diffs, m[:absolute_difference])
    push!(relative_diffs, m[:relative_difference])
    push!(cond_originals, m[:condition_number_original])
    push!(cond_perturbed_vals, m[:condition_number_perturbed])
    push!(max_eigen_diffs, m[:max_eigenvalue_difference])
    push!(max_singular_diffs, m[:max_singular_difference])
    # push!(frob_norms, m[:frobenius_norm])
end

# Plotting the results
f = Figure()
ax = Axis(f[1, 1], xlabel="Noise Level", ylabel="Metric Value", title="Sensitivity Analysis")

plot!(noise_levels, absolute_diffs, label="Absolute Difference", marker=:o)
plot!(noise_levels, relative_diffs, label="Relative Difference", marker=:o)
axislegend(ax, position=:lt)
# plot!(noise_levels, frob_norms, label="Frobenius Norm", marker=:o)
display(f)