如何在Normalized Cuts中用NumPy和Python存储大型矩阵及求特征向量？

一、正确构建并存储大型邻接矩阵 #

100×200图像对应20000个节点，20000×20000的稠密矩阵内存开销极大（约3.2GB），完全不具备可行性，稀疏矩阵是唯一合理的选择，你之前用csr_matrix失败大概率是构建方式有误。

Normalized Cuts的邻接矩阵仅在邻域像素间存在非零值（比如4邻域、8邻域或高斯加权邻域），构建时只需聚焦这些非零元素：

• 先确定邻域规则（比如8邻域包含上下左右及四个对角线像素）
• 遍历每个像素，计算其与邻接像素的相似度（常用高斯相似度：$W_{ij} = e^{{-\frac{(I_i-I_j)}2}{2\sigma^2}}$，$I_i$为像素灰度/颜色特征）
• 用coo_matrix先存储行、列、值的三元组（coo格式最适合增量构建稀疏矩阵），再转换为csr_matrix用于后续计算

示例代码片段：

import numpy as np
from scipy.sparse import coo_matrix, csr_matrix

img = np.random.rand(100, 200)  # 示例输入图像
sigma = 0.1  # 高斯相似度的带宽参数
rows = []
cols = []
data = []

height, width = img.shape
total_nodes = height * width

# 遍历每个像素构建邻接关系
for i in range(height):
    for j in range(width):
        idx = i * width + j
        # 遍历8邻域
        for di in [-1, 0, 1]:
            for dj in [-1, 0, 1]:
                if di == 0 and dj == 0:
                    continue  # 跳过像素自身
                ni, nj = i + di, j + dj
                if 0 <= ni < height and 0 <= nj < width:
                    nidx = ni * width + nj
                    # 计算高斯相似度
                    diff = img[i,j] - img[ni,nj]
                    weight = np.exp(-diff**2 / (2 * sigma**2))
                    rows.append(idx)
                    cols.append(nidx)
                    data.append(weight)

# 构建coo矩阵并转为csr格式
W = coo_matrix((data, (rows, cols)), shape=(total_nodes, total_nodes)).tocsr()
# 计算度矩阵D（对角稀疏矩阵）
D = csr_matrix((np.array(W.sum(axis=1)).flatten(), (range(total_nodes), range(total_nodes))))

二、无需存储完整矩阵求解特征向量的方法 #

Normalized Cuts的核心是求解广义特征值问题：$(D-W)\mathbf{y} = \lambda D\mathbf{y}$，等价于归一化后的标准特征值问题：$L_{sym}\mathbf{z} = \lambda \mathbf{z}$（其中$L_{sym} = D^{-1/2}(D-W)D{-1/2}$，$\mathbf{z}=D^{1/2}\mathbf{y}$）。

不需要显式存储$L_{sym}$或$D-W$，可以用矩阵-向量乘积迭代法（比如Lanczos算法），这类算法仅需每次计算“矩阵×向量”的结果，无需存储完整矩阵：

1. 实现自定义的矩阵-向量乘积函数：将$L_{sym}\mathbf{v}$拆解为三步计算，避免显式构建矩阵
   • 第一步：$\mathbf{u} = D^{-1/2}\mathbf{v}$（元素-wise除以$\sqrt{D_{ii}}$）
   • 第二步：$\mathbf{w} = (D-W)\mathbf{u}$（即$D\mathbf{u} - W\mathbf{u}$，$D\mathbf{u}$为元素乘积，$W\mathbf{u}$可通过稀疏矩阵运算或邻域直接计算）
   • 第三步：$\mathbf{result} = D^{-1/2}\mathbf{w}$
2. 使用scipy.sparse.linalg.eigsh函数，通过matvec参数传入自定义乘积函数，求解最小的几个特征值（Normalized Cuts仅需第二小的特征向量用于分割）

示例代码片段：

from scipy.sparse.linalg import eigsh

# 预计算D的逆平方根（对角元素）
D_inv_sqrt = np.array(1 / np.sqrt(D.diagonal()))

def matvec(v):
    # 计算 D^{-1/2} * v
    u = v * D_inv_sqrt
    # 计算 (D-W)u = Du - Wu
    du = D.dot(u)
    wu = W.dot(u)
    w = du - wu
    # 计算 D^{-1/2} * w
    return w * D_inv_sqrt

# 求解第二小的特征值和特征向量（k=2，取最后一个结果）
vals, vecs = eigsh(matvec, n=total_nodes, k=2, which='SM')
# 转换回原问题的特征向量y
y = vecs[:,1] * D_inv_sqrt

这种方法仅需存储向量和邻域计算的临时变量，内存开销极低，适合处理大尺寸图像。

内容的提问来源于stack exchange，提问作者zo Aoo