Scaling and Equilibration

Scaling and Equilibration

Comprehensive guide for matrix scaling to improve conditioning and solver performance.

When Scaling is Needed

Indicators of Poor Scaling

Matrix Value Ranges

Scaling Methods

Row Scaling

Multiply each row by a scalar: D_r × A

For each row i:
    scale_i = 1 / max_j |A_ij|  (max scaling)
    or
    scale_i = 1 / ||row_i||₂    (norm scaling)
    or
    scale_i = 1 / ||row_i||₁    (sum scaling)

Effect: All row maxima (or norms) become 1.

Preserves: Column space structure.

Column Scaling

Multiply each column by a scalar: A × D_c

For each column j:
    scale_j = 1 / max_i |A_ij|

Effect: All column maxima become 1.

Preserves: Row space structure.

Row and Column Scaling

Apply both: D_r × A × D_c

Balanced approach:

Scaled system: (D_r × A × D_c) × (D_c⁻¹ × x) = D_r × b
Solve for: y = D_c⁻¹ × x
Recover: x = D_c × y

Symmetric Scaling

For symmetric matrices, use same scaling for rows and columns: D × A × D

Preserves symmetry: (DAD)ᵀ = DAD
Essential for CG, MINRES

Computing symmetric scaling:

d_i = 1 / sqrt(|A_ii|)  (diagonal scaling)
or
d_i = 1 / sqrt(max_j |A_ij|)  (max scaling)

Equilibration Algorithms

Ruiz Equilibration

Iterative algorithm to make row and column norms approximately equal:

Repeat until converged:
    For each row i: r_i = ||row_i||_∞
    For each col j: c_j = ||col_j||_∞
    D_r = diag(1/sqrt(r_i))
    D_c = diag(1/sqrt(c_j))
    A = D_r × A × D_c

Properties:

• Converges to doubly stochastic-like scaling
• Usually 5-10 iterations sufficient
• Works for nonsymmetric matrices

Sinkhorn-Knopp (Doubly Stochastic)

For nonnegative matrices, scale to doubly stochastic:

All row sums = 1
All column sums = 1

Algorithm:

Repeat:
    Normalize rows to sum to 1
    Normalize columns to sum to 1
Until converged

Geometric Mean Scaling

Scale by geometric mean of max and min absolute values:

For row i:
    r_max = max_j |A_ij|
    r_min = min_{j: A_ij≠0} |A_ij|
    scale_i = 1 / sqrt(r_max × r_min)

Good for: Matrices with entries spanning many orders of magnitude.

Special Considerations

Preserving Symmetry

CRITICAL: For symmetric matrices, use symmetric scaling only!

Wrong: D_r × A × D_c with D_r ≠ D_c (destroys symmetry)
Right: D × A × D with same D for rows and columns

Check after scaling:

||A_scaled - A_scaled^T||_F / ||A_scaled||_F < ε

Preserving Positive Definiteness

Symmetric scaling preserves definiteness:

If A is SPD and D is nonsingular diagonal:
Then D × A × D is also SPD

Check after scaling (if needed):

Try Cholesky factorization
Or check smallest eigenvalue > 0

Scaling the RHS

When scaling Ax = b → (D_r A D_c) y = D_r b:

Solve: (D_r A D_c) y = D_r b
Recover: x = D_c y
Residual: r = b - Ax = b - A(D_c y)

Important: Scale b with the same D_r used for rows!

Practical Scaling Strategies

Safe Default Strategy

def scale_matrix(A, symmetric=False):
    # Row scaling
    row_norms = np.max(np.abs(A), axis=1)
    D_r = 1.0 / np.maximum(row_norms, 1e-15)

    if symmetric:
        # Use same scaling for rows and columns
        D_c = D_r
    else:
        # Scale by rows first, then columns
        A_scaled = np.diag(D_r) @ A
        col_norms = np.max(np.abs(A_scaled), axis=0)
        D_c = 1.0 / np.maximum(col_norms, 1e-15)

    return D_r, D_c

When to Apply Scaling

Scaling Order

1. Analyze original matrix (symmetry, definiteness)
2. Choose appropriate scaling method
3. Scale matrix and RHS
4. Build preconditioner on scaled system
5. Solve scaled system
6. Unscale solution

Diagnosing Scaling Issues

Before Scaling

Check: ratio of max to min nonzero absolute values
       ||A||_∞ / min_{ij: A_ij≠0} |A_ij|

If ratio > 10⁶: scaling strongly recommended

After Scaling

Verify:
- Row norms approximately equal
- Column norms approximately equal
- No overflow/underflow
- Symmetry preserved (if was symmetric)

Monitoring Convergence

Common Mistakes

Mistake 1: Destroying Symmetry

Problem: Used different row and column scaling on symmetric matrix
Result: CG fails (requires symmetric)
Fix: Use symmetric scaling D × A × D

Mistake 2: Forgetting to Unscale

Problem: Returned y instead of x = D_c × y
Result: Wrong solution
Fix: Always unscale: x = D_c × y

Mistake 3: Not Scaling RHS

Problem: Scaled A but not b
Result: Solving wrong system
Fix: b_scaled = D_r × b

Mistake 4: Scaling Zero Rows

Problem: Row of zeros gets 1/0 = inf scaling
Result: NaN in solution
Fix: Handle zero rows specially (or they indicate singular matrix)

Advanced Topics

Iterative Refinement with Scaling

For very ill-conditioned systems:

1. Scale: Ã = D_r A D_c, b̃ = D_r b
2. Solve: Ã ỹ = b̃ (in lower precision if available)
3. Compute residual: r = b - A(D_c ỹ) (in high precision)
4. Solve: Ã δỹ = D_r r
5. Update: ỹ = ỹ + δỹ
6. Repeat until converged
7. Unscale: x = D_c ỹ

Scaling for Eigenvalue Problems

When solving Ax = λx after scaling D_r A D_c:

Eigenvalues: Same (scaling is similarity transform)
Eigenvectors: v_original = D_c × v_scaled

Block Scaling

For block-structured matrices:

Scale each block independently, or
Scale entire rows/columns of blocks

Useful for: Multiphysics with different variable scales.

Implementation Example

def equilibrate_matrix(A, tol=1e-6, max_iter=20):
    """Ruiz equilibration for general matrices."""
    n, m = A.shape
    D_r = np.ones(n)
    D_c = np.ones(m)

    for _ in range(max_iter):
        # Row scaling
        row_norms = np.max(np.abs(A), axis=1)
        row_norms = np.maximum(row_norms, 1e-15)
        d_r = 1.0 / np.sqrt(row_norms)
        A = np.diag(d_r) @ A
        D_r *= d_r

        # Column scaling
        col_norms = np.max(np.abs(A), axis=0)
        col_norms = np.maximum(col_norms, 1e-15)
        d_c = 1.0 / np.sqrt(col_norms)
        A = A @ np.diag(d_c)
        D_c *= d_c

        # Check convergence
        if np.max(np.abs(row_norms - 1)) < tol:
            if np.max(np.abs(col_norms - 1)) < tol:
                break

    return A, D_r, D_c

def solve_scaled(A, b, solver, symmetric=False):
    """Solve with scaling and unscaling."""
    A_scaled, D_r, D_c = equilibrate_matrix(A)
    b_scaled = D_r * b

    if symmetric:
        # Ensure symmetric scaling
        D = np.sqrt(D_r * D_c)
        A_scaled = np.diag(D) @ A @ np.diag(D)
        b_scaled = D * b
        y = solver(A_scaled, b_scaled)
        x = D * y
    else:
        y = solver(A_scaled, b_scaled)
        x = D_c * y

    return x

Quick Reference