Перекрытие‑осведомлённая 2D‑корреляция Пирсона без нулевого дополнения в SciPy

Как быстро получить карту 2D‑корреляции Пирсона по реальному перекрытию без дополнения нулями: свёрточные суммы через scipy.signal.correlate2d, код и проверка.

Вычисление двумерной взаимной корреляции, которая на каждом смещении ведёт себя как корреляция Пирсона, часто требуется в анализе изображений и сеточных данных датчиков. Сложность в том, чтобы избежать дополнения нулями и использовать только реальное перекрытие массивов для каждого лага. Наивная реализация с вложенными циклами и вызовом pearsonr для каждого лага работает, но чрезвычайно медленна. Ниже — практический способ сделать это быстро, сохранив точное соответствие коэффициенту Пирсона R на области перекрытия при каждом сдвиге.

Проблема

Нужна 2D‑поверхность корреляции формы (двойной размер входа минус один) по каждому измерению — как у scipy.signal.correlate2d в режиме full. Но вместо дополнения нулями у краёв каждая ячейка вывода должна быть коэффициентом корреляции Пирсона, рассчитанным по данным, которые действительно перекрываются на этом смещении. Для входа размера rows × cols величина перекрытия при row_lag и col_lag равна (rows − |row_lag|) × (cols − |col_lag|). Прямолинейное решение — делать срез для каждого лага и вызывать pearsonr на расплющенной области перекрытия, — но оно плохо масштабируется.

Код, который выглядит верным, но не даёт Пирсона для каждого лага

Заманчивый подход — глобально стандартизовать (z‑score) оба массива, один раз запустить correlate2d и поделить на размер перекрытия. Это быстро, но нормализация выполняется глобально, а не для каждого окна перекрытия, поэтому для ненулевых лагов результат не является коэффициентом корреляции Пирсона.

import numpy as np
from scipy.signal import correlate2d

def norm_std(arr):
    return (arr - np.mean(arr)) / np.std(arr)

def xcorr2d_overlap_scaled(u, v):
    assert u.shape == v.shape
    assert u.ndim == 2

    h, w = u.shape
    dr = np.arange(-h + 1, h)
    dc = np.arange(-w + 1, w)

    ov_h = h - np.abs(dr)
    ov_w = w - np.abs(dc)
    overlap = ov_h.reshape((-1, 1)) * ov_w.reshape((1, -1))

    u0 = norm_std(u)
    v0 = norm_std(v)
    raw = correlate2d(u0, v0)
    return raw / overlap

Так получается поверхность, совпадающая с Pearson R только в центре, где перекрытие охватывает весь массив. По мере удаления от центра Пирсон требует локального среднего и локального стандартного отклонения, рассчитанных по окну перекрытия на каждом лаге. Глобальный z‑score этого не учитывает.

Почему наивная нормализация не работает

Корреляция Пирсона на заданном лаге — это ковариация перекрывающихся значений, делённая на произведение их стандартных отклонений; и числитель, и знаменатель считаются по одной и той же области перекрытия. В числителе — сумма попарных произведений минус произведение сумм, масштабированное числом элементов перекрытия. Знаменатель зависит от локальной дисперсии, а та — от локальных сумм и сумм квадратов. Деление глобально нормированной корреляции на размер перекрытия не восстанавливает эти статистики окна и, следовательно, не совпадает с Pearson R, кроме случая полного перекрытия.

Быстрое решение через свёрточные суммы

Ключевая идея — для каждого смещения получить: число перекрывающихся элементов, сумму каждого массива по области перекрытия, сумму квадратов и сумму попарных произведений. Всё это можно получить через correlate2d, свёртывая входы с массивом из единиц той же формы, что и входы. Имея эти величины, можно сформировать числитель и знаменатель Пирсона для каждого лага — за один векторизованный проход.

import numpy as np
from scipy.signal import correlate2d

def pearson_xcorr2d_fast(x, y):
    assert x.shape == y.shape
    x = x.astype(np.float64)
    y = y.astype(np.float64)

    h, w = x.shape
    box = np.ones_like(x)

    overlap_count = correlate2d(box, box)
    sum_x = correlate2d(x, box)
    sum_y = correlate2d(box, y)
    sum_xy = correlate2d(x, y)
    sum_x2 = correlate2d(x * x, box)
    sum_y2 = correlate2d(box, y * y)

    num = sum_xy - (sum_x * sum_y) / overlap_count
    var_x = sum_x2 - (sum_x ** 2) / overlap_count
    var_y = sum_y2 - (sum_y ** 2) / overlap_count
    den = np.sqrt(var_x * var_y)

    with np.errstate(invalid='ignore', divide='ignore'):
        rmap = num / den
    rmap[np.isnan(rmap)] = 0
    return rmap

Выход имеет форму (2*h − 1, 2*w − 1) — как у correlate2d в режиме full. Каждая ячейка — коэффициент корреляции Пирсона двух входов ровно по области, которая перекрывается на данном лаге. Там, где локальная дисперсия равна нулю хотя бы в одном из перекрытий, значение неопределённо; в коде выше такие элементы заменяются нулями.

Проверка промежуточных сумм

Корректность расчёта сумм с учётом перекрытия легко проверить на любом выбранном лаге. Сумма по x для конкретного смещения равна сумме соответствующего среза области перекрытия.

# Пример проверки для положительного смещения
h, w = x.shape
lag_r, lag_c = 3, 7
r_idx = h - 1 + lag_r
c_idx = w - 1 + lag_c
x_slice = x[lag_r:, lag_c:]
np.allclose(sum_x[r_idx, c_idx], np.sum(x_slice))

Это соответствует интуитивной модели: correlate2d с ядром из единиц подсчитывает и агрегирует ровно те элементы, что входят в перекрытие для каждого смещения.

Почему это важно

Для реальных 2D‑данных проход по всем лагам с пересчётом статистик Пирсона на расплющенных окнах непозволительно медленный. Подход на базе свёрток использует correlate2d, чтобы за один раз получить все нужные суммы для каждого смещения, и даёт поверхность, согласующуюся с медленным эталоном. Он также избегает ловушки глобальной нормализации, искажающей значения у краёв. Учтите: для разреженных входов, где в большинстве областей перекрытия дисперсия практически нулевая, Pearson R будет неопределён или равен нулю; тогда как сырая взаимная корреляция из correlate2d всё ещё может показывать пики выравнивания, потому что не нормируется по дисперсии.

Выводы

Используйте «перекрытие‑осведомлённую» корреляцию Пирсона, когда требуется измерять сходство в статистическом смысле на каждом пространственном смещении. Глобальный z‑score плюс корреляция не заменяют нормализацию по каждому лагу. Приём со свёрткой и массивом из единиц позволяет эффективно получить количества, суммы, суммы квадратов и перекрёстные суммы по всем лагам, а из них — корректную поверхность Pearson R. Для проверки сверяйте точечно с простым циклом for на небольшом входе. Если ваши данные по области перекрытия в основном константны, ожидайте нулевые или неопределённые корреляции в этих позициях — для Пирсона это корректно.