COO格式稀疏矩阵进行Permutation置换操作

要对一个COO 格式（Coordinate Format）的稀疏方阵A ∈ R n × n A \in \mathbb{R}^{n \times n}A∈Rn×n应用对称置换（symmetric permutation），即计算：

A ′ = P A P T A' = P A P^TA′=PAPT

其中 ( P ) 是由给定的一维排列数组perm（长度为 ( n )）所定义的置换矩阵，其作用为：

• P [ i , perm [ i ] ] = 1 P[i, \text{perm}[i]] = 1P[i,perm[i]]=1
• 等价地，对向量x xx，有( P x ) i = x perm [ i ] (Px)_i = x_{\text{perm}[i]}(Px)i​=xperm[i]​

那么对矩阵的对称置换等价于：

• 行重排：行i ii变为原来的行perm [ i ] \text{perm}[i]perm[i]
• 列重排：列j jj变为原来的列perm [ j ] \text{perm}[j]perm[j]

对 COO 格式的三元组(row, col, val)，操作非常直接：

✅ 步骤 1：应用对称置换

对每一个非零元素(r, c, v)：

• 新行号 =inv_perm[r]
• 新列号 =inv_perm[c]

为什么是inv_perm？

因为perm[i] = j表示“新位置 i 来自旧位置 j”，即：

• 新矩阵的第 i 行 = 原矩阵的第perm[i]行
  => 原矩阵的第 r 行 → 出现在新矩阵的inv_perm[r]行

所以要将原(r, c)映射到新位置(inv_perm[r], inv_perm[c])

结论：对称置换P A P T P A P^TPAPT对应于将每个 COO 元素的行/列索引通过逆排列映射。

✅ 步骤 2：构造逆排列

给定perm（0-based），其逆排列inv_perm满足：

inv_perm[perm[i]]=i;

代码：

std::vector<int>inv_perm(n);for(inti=0;i<n;++i){inv_perm[perm[i]]=i;}

✅ 步骤 3：应用置换并排序 COO

COO 格式通常要求按(row, col)字典序排序（便于后续转换为 CSR 或其他格式）。

C++ 实现如下（假设使用 0-based 索引）：

#include<vector>#include<algorithm>#include<tuple>structCOOMatrix{intn;// dimensionstd::vector<int>row,col;// row/col indicesstd::vector<double>val;// values (or use template)};COOMatrixapply_symmetric_permutation(constCOOMatrix&A,conststd::vector<int>&perm){intnnz=A.row.size();intn=A.n;// Step 1: build inverse permutationstd::vector<int>inv_perm(n);for(inti=0;i<n;++i){inv_perm[perm[i]]=i;}// Step 2: apply permutation to each entrystd::vector<int>new_row(nnz),new_col(nnz);std::vector<double>new_val=A.val;// copy valuesfor(intk=0;k<nnz;++k){new_row[k]=inv_perm[A.row[k]];new_col[k]=inv_perm[A.col[k]];}// Step 3: sort by (row, col)std::vector<std::tuple<int,int,double>>entries(nnz);for(intk=0;k<nnz;++k){entries[k]=std::make_tuple(new_row[k],new_col[k],new_val[k]);}std::sort(entries.begin(),entries.end());// Step 4: unpack backCOOMatrix A_perm;A_perm.n=n;A_perm.row.resize(nnz);A_perm.col.resize(nnz);A_perm.val.resize(nnz);for(intk=0;k<nnz;++k){std::tie(A_perm.row[k],A_perm.col[k],A_perm.val[k])=entries[k];}returnA_perm;}

🔁 如何“逆置换”？

如果你已经对A AA应用了对称置换得到A ′ = P A P T A' = P A P^TA′=PAPT，那么恢复原矩阵只需应用逆排列的对称置换：

• 原排列为perm，则逆操作的排列是inv_perm
• 即：A = P^T A' P = apply_symmetric_permutation(A_prime, inv_perm)

或者更简单：对A'再次用perm作为排列调用上述函数？❌ 不行！

正确做法：要撤销P A P T P A P^TPAPT，应使用P T A ′ P = ( P − 1 ) A ′ ( P − 1 ) T P^T A' P = (P^{-1}) A' (P^{-1})^TPTA′P=(P−1)A′(P−1)T，而P − 1 P^{-1}P−1对应的排列就是inv_perm。

所以：

COOMatrix A_restored=apply_symmetric_permutation(A_perm,inv_perm_of_perm);

但注意：inv_perm_of_perm就是原始的perm！因为inv_perm的逆是perm。

因此：

• 正向：A1 = apply(..., perm)
• 逆向：A0 = apply(..., inv_perm)，其中inv_perm由perm构造

✅ 补充说明

• 该方法适用于任意稀疏结构（包括非对称矩阵），但对称置换常用于对称矩阵（如刚度矩阵）以改善数值性质。
• COO 输出已按(row, col)排序，可直接用于Eigen::SparseMatrix::setFromTriplets或其他库。
• 若你使用的是size_t或int64_t索引，请相应调整类型。