Examples

This section provides detailed examples of using alsgls for various applications.

Basic SUR Example

Seemingly Unrelated Regressions (SUR) is a classic application of GLS where we have multiple regression equations with correlated error terms.

from alsgls import simulate_sur, als_gls, XB_from_Blist, mse
import numpy as np

# Simulate SUR data
# N_tr=300 observations, K=50 equations, p=3 regressors per equation, k=4 factors
Xs_train, Y_train, Xs_test, Y_test = simulate_sur(
    N_tr=300, N_te=100, K=50, p=3, k=4, seed=0
)

print(f"Training data: {len(Xs_train)} equations")
print(f"Equation shapes: {[X.shape for X in Xs_train[:3]]}...")
print(f"Response matrix: {Y_train.shape}")

# Fit ALS model
B, F, D, memory_usage, info = als_gls(Xs_train, Y_train, k=4, sweeps=10, rel_tol=1e-6)

print(f"Converged in {len(info['nll_trace']) - 1} iterations")
print(f"Final objective: {info['nll_trace'][-1]:.6f}")
print(f"Memory usage: {memory_usage:.3f} MB")

# Make predictions and evaluate
Y_pred = XB_from_Blist(Xs_test, B)
test_mse = mse(Y_test, Y_pred)
print(f"Test MSE: {test_mse:.6f}")

Large-Scale Example

This example demonstrates the memory advantages for larger problems:

from alsgls import simulate_sur, als_gls
import psutil
import os


def get_memory_usage():
    """Get current memory usage in MB"""
    process = psutil.Process(os.getpid())
    return process.memory_info().rss / 1024 / 1024


# Large problem: 200 equations, 500 observations
print("Generating large-scale data...")
Xs_tr, Y_tr, Xs_te, Y_te = simulate_sur(N_tr=500, N_te=200, K=200, p=4, k=6)

print(f"Problem size: {Y_tr.shape[0]} obs × {Y_tr.shape[1]} equations")
print(f"Total parameters in dense Σ: {Y_tr.shape[1] ** 2:,}")
print(f"Parameters in low-rank model: {Y_tr.shape[1] * 6 + Y_tr.shape[1]:,}")

mem_before = get_memory_usage()

print("\\nFitting ALS model...")
B_als, F_als, D_als, mem_als, info_als = als_gls(Xs_tr, Y_tr, k=6, sweeps=8)
print(f"ALS converged in {len(info_als['nll_trace']) - 1} iterations")


Y_pred_als = XB_from_Blist(Xs_te, B_als)
mse_als = mse(Y_te, Y_pred_als)

dense_mb = Y_tr.shape[1] ** 2 * 8 / 1024 / 1024
print(f"\\nResults:")
print(f"ALS - Memory: {mem_als:.1f}MB, MSE: {mse_als:.6f}")
print(f"A dense Sigma would need: {dense_mb:.1f}MB")

Custom Data Example

Working with your own data instead of simulated data:

import numpy as np
from alsgls import als_gls, XB_from_Blist

# Prepare your data
# Xs should be a list of design matrices, one per equation
# Y should be an (N, K) matrix of responses

# Example: 3 equations with different numbers of regressors
N = 200  # observations
K = 3  # equations

# Generate some example data
np.random.seed(42)

# Equation 1: 2 regressors
X1 = np.random.randn(N, 2)

# Equation 2: 3 regressors
X2 = np.random.randn(N, 3)

# Equation 3: 1 regressor
X3 = np.random.randn(N, 1)

Xs = [X1, X2, X3]

# Correlated responses (you would load your actual data here)
true_B = [
    np.array([[1.5], [-0.8]]),  # coeffs for eq 1
    np.array([[0.5], [1.2], [-0.3]]),  # coeffs for eq 2
    np.array([[2.0]]),
]  # coeffs for eq 3

# Generate responses with correlation structure
Y = np.column_stack(
    [
        X1 @ true_B[0].flatten() + 0.1 * np.random.randn(N),
        X2 @ true_B[1].flatten() + 0.1 * np.random.randn(N),
        X3 @ true_B[2].flatten() + 0.1 * np.random.randn(N),
    ]
)

# Add some cross-equation correlation
factor = np.random.randn(N, 1)
Y += 0.3 * factor @ np.random.randn(1, K)

print("Data shapes:")
for i, X in enumerate(Xs):
    print(f"  Equation {i + 1}: X{i + 1} = {X.shape}")
print(f"  Responses: Y = {Y.shape}")

# Fit the model
B_hat, F, D, memory, info = als_gls(Xs, Y, k=2, sweeps=10)

print(f"\\nEstimated coefficients:")
for i, b in enumerate(B_hat):
    print(f"  Equation {i + 1}: {b.flatten()}")
    print(f"  True values: {true_B[i].flatten()}")
    print(f"  Error: {np.linalg.norm(b.flatten() - true_B[i].flatten()):.4f}")

Factor Structure Analysis

Examining the estimated factor structure:

from alsgls import simulate_sur, als_gls
import matplotlib.pyplot as plt

# Simulate data with known factor structure
Xs_tr, Y_tr, _, _ = simulate_sur(N_tr=400, N_te=100, K=30, p=3, k=3)

# Fit model
B, F, D, _, _ = als_gls(Xs_tr, Y_tr, k=3, sweeps=15)

print(f"Factor loadings shape: {F.shape}")
print(f"Diagonal variances shape: {D.shape}")

# Examine factor loadings
print("\\nFactor loadings (first 5 equations):")
print(F[:5, :])

# Compute explained variance by factors
factor_var = np.var(F @ F.T, axis=1)
total_var = factor_var + D
explained_ratio = factor_var / total_var

print(f"\\nVariance explained by factors:")
print(f"  Mean: {explained_ratio.mean():.3f}")
print(f"  Min:  {explained_ratio.min():.3f}")
print(f"  Max:  {explained_ratio.max():.3f}")

# Plot factor loadings heatmap
plt.figure(figsize=(8, 6))
plt.imshow(F, aspect="auto", cmap="RdBu_r")
plt.colorbar(label="Loading")
plt.xlabel("Factor")
plt.ylabel("Equation")
plt.title("Factor Loadings Matrix")
plt.show()

Cross-Validation Example

Selecting the optimal number of factors using cross-validation:

from alsgls import simulate_sur, als_gls, XB_from_Blist, mse
import matplotlib.pyplot as plt
import numpy as np

# Generate data
Xs_tr, Y_tr, Xs_val, Y_val = simulate_sur(N_tr=300, N_te=150, K=40, p=3, k=4)

# Test different numbers of factors
k_values = range(1, 11)
val_mses = []

for k in k_values:
    print(f"Testing k={k}...")

    # Fit model
    B, F, D, _, info = als_gls(Xs_tr, Y_tr, k=k, sweeps=10)

    # Validate
    Y_pred = XB_from_Blist(Xs_val, B)
    val_mse = mse(Y_val, Y_pred)
    val_mses.append(val_mse)

    print(f"  Validation MSE: {val_mse:.6f}")

# Find optimal k
best_k = k_values[np.argmin(val_mses)]
print(f"\\nOptimal number of factors: k={best_k}")

# Plot validation curve
plt.figure(figsize=(8, 5))
plt.plot(k_values, val_mses, "o-")
plt.axvline(best_k, color="red", linestyle="--", label=f"Optimal k={best_k}")
plt.xlabel("Number of factors (k)")
plt.ylabel("Validation MSE")
plt.title("Factor Selection via Cross-Validation")
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

Performance Profiling

Detailed timing and memory profiling:

import time
from alsgls import simulate_sur, als_gls


def profile_solver(solver_func, Xs, Y, k, **kwargs):
    """Profile a solver function"""
    start_time = time.time()
    result = solver_func(Xs, Y, k=k, **kwargs)
    end_time = time.time()

    B, F, D, memory, info = result
    runtime = end_time - start_time

    return {
        "runtime": runtime,
        "memory": memory,
        "iterations": len(info["nll_trace"]) - 1,
        "objective": info["nll_trace"][-1],
        "B": B,
        "F": F,
        "D": D,
    }


# Test different problem sizes
problem_sizes = [(100, 30), (200, 60), (300, 90)]

for N, K in problem_sizes:
    print(f"\\nProblem size: N={N}, K={K}")

    # Generate data
    Xs, Y, _, _ = simulate_sur(N_tr=N, N_te=50, K=K, p=3, k=5)

    # Profile ALS
    als_result = profile_solver(als_gls, Xs, Y, k=5, sweeps=8)

    print(
        f"ALS: {als_result['runtime']:.2f}s, {als_result['memory']:.1f}MB, "
        f"{als_result['iterations']} iter"
    )
    print(f"Dense Sigma would need: {K**2 * 8 / 1024 / 1024:.1f}MB")