我有一个 NxN 非负
Eigen::MatrixXicost_matrixEigen::VectorXirowVectorcolVectorcost_matrix(i, j)rowVector(i)colVector(j)我知道有一些简单的解决方案，例如遍历所有元素，但它们花费了太多时间。我想要 Eigen C++ 中最有效的方法。
这是我当前的代码，它比使用
Eigen::DenseBase::visitEigen::MatrixXi cost_matrix(4,4);
cost_matrix<<1, 2, 3, 0,
5, 0, 0, 8,
9, 8, 7, 6,
0, 2, 1, 5;
Eigen::VectorXi rowCover(4);
Eigen::VectorXi colCover(4);
rowCover << 0, 0, 1, 1;
colCover << 1, 1, 0, 1;
//A data demo. Cost_matrix is a NxN nonnegative int matrix.
//RowCover and colCover are N nonnegative int vector.
int i, j;
cost_matrix.colwise() += rowCover;
cost_matrix.rowwise() += colCover.transpose();
Eigen::Index zeroRow, zeroCol;
int min = find_zero_matrix.transpose().minCoeff(&zeroCol, &zeroRow);
i = zeroRow;
j = zeroCol;
//the result should be
//i = 1
//j = 2
正如 @chtz 已经在评论中建议的那样，Eigen 在这里不会真正帮助你。我们仍然可以尝试找到更快的版本。
这是我的想法：
首先，我扩展/修复您的代码以获得工作参考实现。
/**
* Reference implementation
*
* Finds lowest row and column (prioritizes low row number)
* where cost_matrix is zero and both row_cover and col_cover are zero
*
* Returns {-1, -1} if no result is found.
* The input arrays can be overwritten by the implementation as desired
*/
std::pair<int, int> reference(
Eigen::MatrixXi& cost_matrix,
Eigen::VectorXi& row_cover,
Eigen::VectorXi& col_cover) noexcept
{
int i, j;
cost_matrix.colwise() += row_cover;
cost_matrix.rowwise() += col_cover.transpose();
Eigen::Index zeroRow, zeroCol;
int min = cost_matrix.transpose().minCoeff(&zeroCol, &zeroRow);
i = zeroRow;
j = zeroCol;
if(min) // no result found
return {-1, -1};
return {i, j};
}
遍历整个输入矩阵似乎非常浪费，除非覆盖向量通常为零而输入矩阵不是。如果封面包含许多非零条目，最好将它们转换为零索引的向量。
/**
* Overwrite the front of cover with indices where it was zero
*
* The return value gives the total number of zeros found.
* The values behind that are left in undefined state
*/
static Eigen::Index to_zero_indices(Eigen::VectorXi& cover) noexcept
{
Eigen::Index outpos = 0;
for(Eigen::Index i = 0, n = cover.size(); i < n; ++i) {
/* Loop body is written to be branchless */
int value = cover[i];
cover[outpos] = static_cast<int>(i);
outpos += ! value;
}
return outpos;
}
现在我们只需要检查矩阵中在两个cover中都为零的条目。
/**
* Alternative implementation based on sparse probing of the cost matrix
*
* API-compatible to reference but in practice changes the cover vectors,
* not the cost matrix
*/
std::pair<int, int> probe_indices(
Eigen::MatrixXi& cost_matrix,
Eigen::VectorXi& row_cover,
Eigen::VectorXi& col_cover) noexcept
{
const Eigen::Index rows_n = to_zero_indices(row_cover);
const Eigen::Index cols_n = to_zero_indices(col_cover);
for(int i: row_cover.head(rows_n))
for(int j: col_cover.head(cols_n))
if(! cost_matrix(i, j))
return {i, j};
return {-1, -1};
}
为了获得最佳结果，请使用
-DNDEBUG矩阵大小 N 约为 20，但每帧会运行很多次。 0 的百分比约为 1/N（可能每行有几个 0）
我假设这个 1/N 不适用于覆盖向量。否则，即使找到一个条目，您的机会也很小。我随意决定用 4/N 来测试封面，用 1/N 来测试成本矩阵。
完整的测试和基准测试结果如下。在我的系统上，使用所选参数集，我的版本快了大约 10 倍，即使我将覆盖率更改为全零，速度也快了 7.5 倍。即使在绝对最坏的情况下——覆盖全为零，矩阵全为——它的速度仍然是原来的两倍。
#include <Eigen/Dense>
#include <chrono>
#include <iostream>
#include <random>
// using std::default_random_engine, std::bernoulli_distribution
#include <utility>
// using std::pair
/**
* Fill the output with zeros and ones depending on the percentage
*
* Results are random but reproducible across process executions.
* Input percentage is expected to be between 0 and 1
*/
void fill_pct(Eigen::Ref<Eigen::MatrixXi> out, double pct_ones)
{
static std::default_random_engine::result_type seed = 0xdeadbeef;
std::default_random_engine rng {seed++};
std::bernoulli_distribution distr {pct_ones};
out = Eigen::MatrixXi::NullaryExpr(
out.rows(), out.cols(),
[&]() -> int { return distr(rng); });
}
int main()
{
using clock_t = std::chrono::steady_clock;
const int size = 20;
const double pct_zero = 1. / size;
const double pct_zero_cover = 4. / size;
const int repetitions = 1000000;
Eigen::MatrixXi cost_matrix, matrix_copy;
Eigen::VectorXi row_cover, col_cover, row_copy, col_copy;
clock_t::duration reference_time {}, bench_time {};
for(int i = 0; i < repetitions; ++i) {
cost_matrix.resize(size, size);
fill_pct(cost_matrix, 1. - pct_zero);
matrix_copy = cost_matrix;
for(Eigen::VectorXi* cover: {&row_cover, &col_cover}) {
cover->resize(size);
fill_pct(*cover, 1. - pct_zero_cover);
}
row_copy = row_cover;
col_copy = col_cover;
clock_t::time_point t1 = clock_t::now();
std::pair<int, int> bench_result = probe_indices(
cost_matrix, row_cover, col_cover);
clock_t::time_point t2 = clock_t::now();
std::pair<int, int> ref_result = reference(
matrix_copy, row_copy, col_copy);
clock_t::time_point t3 = clock_t::now();
bench_time += t2 - t1;
reference_time += t3 - t2;
if(bench_result != ref_result)
std::cout << bench_result.first << " != " << ref_result.first
<< ' ' << bench_result.second << " != " << ref_result.second
<< '\n';
}
using std::chrono::milliseconds;
std::cout << "Reference " << std::chrono::duration_cast<milliseconds>(
reference_time).count() * 1e-3
<< "\nBench " << std::chrono::duration_cast<milliseconds>(
bench_time).count() * 1e-3
<< '\n';
}
附录： 我使用显式 AVX512 矢量化进行了一些测试。首先让
to_zero_indicesvpcompressd