Linear Algebra and Numerical Computing Questions

When summing many small floating-point numbers you can lose precision because rounding errors accumulate. Useful approaches: pairwise (divide-and-conquer) summation, compensated summation (Kahan or Neumaier), and using higher precision (float128/longdouble or Decimal). Pairwise reduces error growth; Kahan/Neumaier corrects running error; higher precision reduces rounding but costs speed/memory.

Implementations and trade-offs:

1. Pairwise (recursive, numerically stable, vectorizable)

• O(n) time, log n recursion depth, good accuracy, works well with NumPy arrays.

import numpy as np

def pairwise_sum(x):
    x = np.asarray(x, dtype=np.float64)
    n = x.size
    if n == 0:
        return 0.0
    if n == 1:
        return float(x[0])
    mid = n // 2
    return pairwise_sum(x[:mid]) + pairwise_sum(x[mid:])

2. Kahan / Neumaier (online, O(n), Python loop — slower but improved accuracy)

def kahan_sum(x):
    s = 0.0
    c = 0.0
    for v in x:
        y = float(v) - c
        t = s + y
        c = (t - s) - y
        s = t
    return s

def neumaier_sum(x):
    s = 0.0
    c = 0.0
    for v in x:
        t = s + v
        if abs(s) >= abs(v):
            c += (s - t) + v
        else:
            c += (v - t) + s
        s = t
    return s + c

3. Higher precision

• Use longdouble when available: np.sum(arr, dtype=np.longdouble) then cast back. Or Python's decimal for exactness (much slower).

Group-sum in pandas / NumPy:

• For many groups, vectorize by aggregating arrays per group then using pairwise or Neumaier per group. Example using pandas groupby + apply:

import pandas as pd

def stable_group_sum(df, group_col, value_col, method='pairwise'):
    grp = df.groupby(group_col)[value_col]
    if method == 'pairwise':
        return grp.apply(lambda s: pairwise_sum(s.values))
    elif method == 'neumaier':
        return grp.apply(lambda s: neumaier_sum(s.values))
    elif method == 'numpy_longdouble':
        return grp.apply(lambda s: float(np.sum(s.values, dtype=np.longdouble)))
    else:
        raise ValueError("unknown method")

Practical guidance:

• If performance matters and arrays are large, prefer pairwise (recursive) implemented in C (NumPy’s sum often uses pairwise/tree reduction already; verify with your NumPy build).
• Use Kahan/Neumaier for streaming sums or when you can’t rely on library reductions.
• Use longdouble or compensated sums when final accuracy is critical.
• Benchmark: Python loops (Kahan) are slower; prefer vectorized/pure-C solutions for large datasets.

NumPy broadcasting lets arrays with different shapes participate in arithmetic by virtually expanding the smaller array along dimensions of size 1 (without copying data) so shapes become compatible. Rules: align trailing dimensions, dimensions must be equal or one of them is 1, or the dimension is missing (treated as 1).

Concrete example: add a (3,1) column vector to a (3,4) matrix — the (3,1) is broadcast across the 4 columns.

import numpy as np

A = np.array([[1],[2],[3]])        # shape (3,1)
B = np.array([[10,20,30,40],
              [11,21,31,41],
              [12,22,32,42]])       # shape (3,4)

C = A + B                           # A is broadcast to shape (3,4)
print(C)
# Equivalent using Python loop (no broadcasting):
C_loop = np.empty_like(B)
for i in range(A.shape[0]):
    for j in range(B.shape[1]):
        C_loop[i,j] = A[i,0] + B[i,j]

Why broadcasting is faster for ML preprocessing:

• Uses vectorized C loops in NumPy, avoiding Python-level loops and interpreter overhead.
• No extra memory copy for repeating values (views), so lower memory bandwidth.
• Enables concise, readable code and lets BLAS/NumPy optimize memory access — critical when applying transforms to large datasets during feature scaling, batching, or augmentation.

To detect outliers using Median Absolute Deviation (MAD) in a fully vectorized way, compute the global median ignoring NaNs, the MAD (median of absolute deviations, ignoring NaNs), convert differences to a "modified z‑score" using the 0.6745 constant (so for normal data MAD estimates the std), and mark points with |z| > thresh as outliers. Handle MAD==0 and all-NaN arrays safely and return a boolean mask same shape as input.

import numpy as np

def mad_outlier_mask(x: np.ndarray, thresh: float = 3.5) -> np.ndarray:
    """
    Return boolean mask of same shape as x where True indicates an outlier
    according to the modified z-score based on MAD. NaNs in x are ignored
    when computing medians; NaN positions in output are False.
    """
    x = np.asarray(x)
    # Flatten for global stats, but remember shape
    flat = x.ravel()
    # Compute median ignoring NaNs
    med = np.nanmedian(flat)
    # Absolute deviations from median (NaNs remain NaN)
    abs_dev = np.abs(flat - med)
    # MAD: median of absolute deviations, ignoring NaNs
    mad = np.nanmedian(abs_dev)

    # If MAD is zero or nan (e.g., all values equal or all-NaN), no robust scaling possible.
    # In that case return all-False mask (no outliers) except keep NaN positions as False.
    if not np.isfinite(mad) or mad == 0:
        return np.zeros_like(x, dtype=bool)

    # Modified z-score: 0.6745 scales MAD to be comparable to std for normal data
    modified_z = 0.6745 * (flat - med) / mad
    mask_flat = np.abs(modified_z) > thresh
    # Ensure NaN positions are False (not marked as outlier)
    mask_flat = np.where(np.isnan(flat), False, mask_flat)

    return mask_flat.reshape(x.shape)

Key points:

• Vectorized using numpy only; no Python loops.
• Ignores NaNs in medians via np.nanmedian; NaN positions in mask are False.
• Uses modified z-score = 0.6745*(x - median)/MAD and threshold default 3.5 (common convention).
  Time complexity: O(n) where n = number of elements. Space: O(n) extra for flattened arrays.
  Edge cases:
• All-NaN input -> returns all-False mask.
• MAD == 0 (identical values) -> returns all-False (no outliers).
  Alternatives:
• Compute MAD per-axis for multidimensional arrays.
• If MAD==0 but variance exists, fallback to std-based z-score.

Approach: use NumPy broadcasting to compute (a_i - b_j)^2 without Python loops.

import numpy as np

def sq_euclidean_matrix(a, b):
    """
    a: shape (n,)
    b: shape (m,)
    returns: D shape (n, m) where D[i,j] = (a[i] - b[j])**2
    """
    a = np.asarray(a)
    b = np.asarray(b)
    # broadcasting: a[:, None] has shape (n,1), b[None, :] has shape (1,m)
    return (a[:, None] - b[None, :])**2

Key points:

• This computes the full (n, m) matrix in one expression using vectorized operations (no Python loops).
• Time complexity: O(nm). Space: O(nm) for the output plus transient memory for the broadcast expression (NumPy computes the difference into a new (n,m) array then squares into it or a separate temp depending on implementation).

Memory cost and reduction technique:

• If n*m is large, the full (n,m) array may not fit in RAM. Two ways to reduce memory:
  1. Compute in chunks along one axis (generator or loop over blocks) and write results to disk / downstream consumer:
     python

     def sq_euclidean_matrix_chunked(a, b, block=1000):
         n = a.shape[0]
         for i in range(0, n, block):
             yield (a[i:i+block, None] - b[None, :])**2  # yields blocks of shape (block, m)
     

  2. Use the algebraic identity to avoid a large temporary: (a - b)^2 = a^2 + b^2 - 2ab. Precompute a_sq and b_sq then compute result (still O(n*m) output but can combine with chunking to reduce peak memory).
• Also consider using lower-precision dtype (float32) or memory-mapped arrays (np.memmap) for very large matrices.