我有一个矩阵NxM(通常为10k X 10k元素),描述了一个基础集合。每一行代表一个对象,每一列代表一个特定的特征。例如,在这个矩阵中:
f1 f2 f3
x1 0 4 -1
x2 1 0 5
x3 4 0 0
x4 0 1 0
对象 x1 在特征 1 中的值为 0,在特征 2 中的值为 4,在特征 -1 中的值为 0。这些值是一般实数(双精度)。
我需要计算所有对象(所有行对)之间的多个自定义距离/不相似度。为了比较,我想计算 L1(曼哈顿)和L2(欧几里得)距离。
我使用 Eigen 库执行大部分计算。为了计算 L2(欧几里得),我使用以下观察结果:对于大小为 n 的两个向量 a 和 b,我们有:
||a - b||^2 = (a_1 - b_1)^2 + (a_2 - b_2)^2 + ... +(a_n - b_n)^2 = a_1^2 + b_1^2 - 2 a_1 b_1 + a_2^2 + b_2^2 - 2 a_2 b_2 + ... + a_n^2 + b_n^2 - 2 a_n b_n = a . a + b . b - 2ab
换句话说,我们使用向量自身的点积重新编写平方范数,并减去它们之间的两倍点积。从那里,我们只需取平方即可完成。有一次,我很久以前发现了这个技巧,不幸的是我丢失了作者的参考。
无论如何,这使得我们可以使用 Eigen(在C++中)编写一段漂亮的代码:
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> matrix, XX, D;
// Load matrix here, for example
// matrix << 0, 4, -1,
// 1, 0, 5,
// 4, 0, 0,
// 0, 1, 0;
const auto N = matrix.rows();
XX.resize(N, 1);
D.resize(N, N);
XX = matrix.array().square().rowwise().sum();
D.noalias() = XX * Eigen::MatrixXd::Ones(1, N) +
Eigen::MatrixXd::Ones(N, 1) * XX.transpose();
D -= 2 * matrix * matrix.transpose();
D = D.cwiseSqrt();
对于10k X 10k的矩阵,我们能够在不到1分钟的时间内计算出所有对象/行之间的L2距离(使用2个核心/4个线程),我认为这对于我的目的来说是良好的性能。Eigen能够组合操作并使用多个低/高级优化来执行这些计算。在这种情况下,Eigen使用了所有可用的核心(当然,我们可以进行调整)。
然而,我仍然需要计算L1距离,但我无法想出一个与Eigen配合使用且获得良好性能的代数形式。到目前为止,我只有以下内容:
const auto N = matrix.rows();
for(long i = 0; i < N - 1; ++i) {
const auto &row = matrix.row(i);
#ifdef _OPENMP
#pragma omp parallel for shared(row)
#endif
for(long j = i + 1; j < N; ++j) {
distance(i, j) = (row - matrix.row(j)).lpNorm<1>();
}
}
但是这需要更长的时间:对于相同的10k X 10k矩阵,此代码使用3.5分钟,考虑到L1和L2在其原始形式中非常接近,这要差得多:
L1(a, b) = sum_i |a_i - b_i|
L2(a, b) = sqrt(sum_i |a_i - b_i|^2)
有没有什么想法可以如何改变L1,使用Eigen进行快速计算?或者有更好的方法来做到这一点,我只是没有想出来。
非常感谢您的帮助!
卡洛斯
D = XX.replicate(1,n) + XX.transpose().replicate(n,1);
- ggael