Source

`tests/benchmarks/eigen_compare_bench.cpp`

2// Original work; see ACKNOWLEDGMENTS.md for the open-source ideas we build upon.

3// cheatah::linalg vs Eigen — a native, single-threaded C++ dense-linear-algebra

4// face-off (Eigen is the reference C++ library). Both sides run the SAME operation on

5// the SAME deterministic data at the SAME sizes the NumPy comparison uses, built at

6// `-O3 -march=native`. Eigen is held to one thread (it is single-threaded by default

7// without OpenMP anyway) so this is an apples-to-apples per-core comparison — the

8// operating point cheatah is designed for.

9//

10// Pairs are named BM_<op>_cheatah / BM_<op>_eigen so the two rows sit together in the

11// report. Build: cmake --preset release-benchmarks (CHEATAH_BUILD_BENCHMARKS=ON).

12#include "linalg.hpp"

13#include "ndarray.hpp"

15#include <cstdint>

16#include <vector>

18#include <Eigen/Dense>

19#include <benchmark/benchmark.h>

21namespace nd = cheatah::ndarray;

22namespace la = cheatah::linalg;

24// Hold Eigen to a single thread so this is a strict per-core comparison (cheatah's

25// numeric core is single-threaded by design). Without OpenMP Eigen is single-threaded

26// anyway — this just makes it explicit and defensive. Runs before main().

27static const int kEigenThreads = (Eigen::setNbThreads(1), Eigen::nbThreads());

29namespace {

31// Deterministic, diagonally-dominant (well-conditioned) matrix entry.

32double gen(long long i, long long j, long long n) {

33 double v = static_cast<double>(((i * 131 + j * 97 + 1) % 7) - 3);

34 if (i == j) v += static_cast<double>(10 * n);

35 return v;

36}

38// Build the same n×n matrix as a cheatah NDArray…

39nd::NDArray cheatah_matrix(long long n) {

40 std::vector<double> buf(static_cast<std::size_t>(n * n));

41 for (long long i = 0; i < n; ++i)

42 for (long long j = 0; j < n; ++j) buf[static_cast<std::size_t>(i * n + j)] = gen(i, j, n);

43 return nd::reshape(nd::array(buf), {n, n});

44}

45// …and as an Eigen matrix.

46Eigen::MatrixXd eigen_matrix(long long n) {

47 Eigen::MatrixXd m(n, n);

48 for (long long i = 0; i < n; ++i)

49 for (long long j = 0; j < n; ++j) m(i, j) = gen(i, j, n);

50 return m;

51}

53// A symmetric well-conditioned matrix (for the symmetric eigensolver).

54double sym(long long i, long long j, long long n) {

55 const long long a = i < j ? i : j, b = i < j ? j : i;

56 double v = static_cast<double>(((a * 53 + b * 29 + 1) % 7) - 3);

57 if (i == j) v += static_cast<double>(10 * n);

58 return v;

59}

60nd::NDArray cheatah_sym(long long n) {

61 std::vector<double> buf(static_cast<std::size_t>(n * n));

62 for (long long i = 0; i < n; ++i)

63 for (long long j = 0; j < n; ++j) buf[static_cast<std::size_t>(i * n + j)] = sym(i, j, n);

64 return nd::reshape(nd::array(buf), {n, n});

65}

66Eigen::MatrixXd eigen_sym(long long n) {

67 Eigen::MatrixXd m(n, n);

68 for (long long i = 0; i < n; ++i)

69 for (long long j = 0; j < n; ++j) m(i, j) = sym(i, j, n);

70 return m;

71}

73// ---- dot (length sweep) -------------------------------------------------

74void BM_dot_cheatah(benchmark::State& state) {

75 const long long n = state.range(0);

76 const nd::NDArray a = nd::full({n}, 1.0), b = nd::full({n}, 2.0);

77 for (auto _ : state) benchmark::DoNotOptimize(la::dot(a, b));

78}

79void BM_dot_eigen(benchmark::State& state) {

80 const long long n = state.range(0);

81 const Eigen::VectorXd a = Eigen::VectorXd::Constant(n, 1.0), b = Eigen::VectorXd::Constant(n, 2.0);

82 for (auto _ : state) benchmark::DoNotOptimize(a.dot(b));

83}

84BENCHMARK(BM_dot_cheatah)->Arg(64)->Arg(16384);

85BENCHMARK(BM_dot_eigen)->Arg(64)->Arg(16384);

87// ---- matmul -------------------------------------------------------------

88void BM_matmul_cheatah(benchmark::State& state) {

89 const long long n = state.range(0);

90 const nd::NDArray a = cheatah_matrix(n), b = cheatah_matrix(n);

91 for (auto _ : state) benchmark::DoNotOptimize(la::matmul(a, b));

92}

93void BM_matmul_eigen(benchmark::State& state) {

94 const long long n = state.range(0);

95 const Eigen::MatrixXd a = eigen_matrix(n), b = eigen_matrix(n);

96 for (auto _ : state) { Eigen::MatrixXd c = a * b; benchmark::DoNotOptimize(c.data()); }

97}

98BENCHMARK(BM_matmul_cheatah)->Arg(32)->Arg(96);

99BENCHMARK(BM_matmul_eigen)->Arg(32)->Arg(96);

100

101// ---- solve (A x = b) ----------------------------------------------------

102void BM_solve_cheatah(benchmark::State& state) {

103 const long long n = state.range(0);

104 const nd::NDArray a = cheatah_matrix(n), b = nd::full({n}, 1.0);

105 for (auto _ : state) benchmark::DoNotOptimize(la::solve(a, b));

106}

107void BM_solve_eigen(benchmark::State& state) {

108 const long long n = state.range(0);

109 const Eigen::MatrixXd a = eigen_matrix(n);

110 const Eigen::VectorXd b = Eigen::VectorXd::Constant(n, 1.0);

111 for (auto _ : state) { Eigen::VectorXd x = a.partialPivLu().solve(b); benchmark::DoNotOptimize(x.data()); }

112}

113BENCHMARK(BM_solve_cheatah)->Arg(32)->Arg(64);

114BENCHMARK(BM_solve_eigen)->Arg(32)->Arg(64);

115

116// ---- inv ----------------------------------------------------------------

117void BM_inv_cheatah(benchmark::State& state) {

118 const long long n = state.range(0);

119 const nd::NDArray a = cheatah_matrix(n);

120 for (auto _ : state) benchmark::DoNotOptimize(la::inv(a));

121}

122void BM_inv_eigen(benchmark::State& state) {

123 const long long n = state.range(0);

124 const Eigen::MatrixXd a = eigen_matrix(n);

125 for (auto _ : state) { Eigen::MatrixXd c = a.inverse(); benchmark::DoNotOptimize(c.data()); }

126}

127BENCHMARK(BM_inv_cheatah)->Arg(32)->Arg(64);

128BENCHMARK(BM_inv_eigen)->Arg(32)->Arg(64);

129

130// ---- det ----------------------------------------------------------------

131void BM_det_cheatah(benchmark::State& state) {

132 const long long n = state.range(0);

133 const nd::NDArray a = cheatah_matrix(n);

134 for (auto _ : state) benchmark::DoNotOptimize(la::det(a));

135}

136void BM_det_eigen(benchmark::State& state) {

137 const long long n = state.range(0);

138 const Eigen::MatrixXd a = eigen_matrix(n);

139 for (auto _ : state) benchmark::DoNotOptimize(a.determinant());

140}

141BENCHMARK(BM_det_cheatah)->Arg(64);

142BENCHMARK(BM_det_eigen)->Arg(64);

143

144// ---- eigvalsh (symmetric eigenvalues) -----------------------------------

145void BM_eigvalsh_cheatah(benchmark::State& state) {

146 const long long n = state.range(0);

147 const nd::NDArray a = cheatah_sym(n);

148 for (auto _ : state) benchmark::DoNotOptimize(la::eigvalsh(a));

149}

150void BM_eigvalsh_eigen(benchmark::State& state) {

151 const long long n = state.range(0);

152 const Eigen::MatrixXd a = eigen_sym(n);

153 for (auto _ : state) {

154 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(a, Eigen::EigenvaluesOnly);

155 benchmark::DoNotOptimize(es.eigenvalues().data());

156 }

157}

158BENCHMARK(BM_eigvalsh_cheatah)->Arg(8)->Arg(64);

159BENCHMARK(BM_eigvalsh_eigen)->Arg(8)->Arg(64);

160

161// ---- svdvals (singular values only) -------------------------------------

162void BM_svdvals_cheatah(benchmark::State& state) {

163 const long long n = state.range(0);

164 const nd::NDArray a = cheatah_matrix(n);

165 for (auto _ : state) benchmark::DoNotOptimize(la::svdvals(a));

166}

167void BM_svdvals_eigen(benchmark::State& state) {

168 const long long n = state.range(0);

169 const Eigen::MatrixXd a = eigen_matrix(n);

170 for (auto _ : state) {

171 Eigen::BDCSVD<Eigen::MatrixXd> svd(a); // no U/V computed -> singular values only

172 benchmark::DoNotOptimize(svd.singularValues().data());

173 }

174}

175BENCHMARK(BM_svdvals_cheatah)->Arg(64);

176BENCHMARK(BM_svdvals_eigen)->Arg(64);

177

178// ---- svd (full U + singular values + Vᵀ) --------------------------------

179void BM_svd_cheatah(benchmark::State& state) {

180 const long long n = state.range(0);

181 const nd::NDArray a = cheatah_matrix(n);

182 for (auto _ : state) benchmark::DoNotOptimize(la::svd(a).s);

183}

184void BM_svd_eigen(benchmark::State& state) {

185 const long long n = state.range(0);

186 const Eigen::MatrixXd a = eigen_matrix(n);

187 for (auto _ : state) {

188 Eigen::BDCSVD<Eigen::MatrixXd> svd(a, Eigen::ComputeFullU | Eigen::ComputeFullV);

189 benchmark::DoNotOptimize(svd.singularValues().data());

190 }

191}

192BENCHMARK(BM_svd_cheatah)->Arg(64)->Arg(96);

193BENCHMARK(BM_svd_eigen)->Arg(64)->Arg(96);

194

195// ---- cholesky -----------------------------------------------------------

196void BM_cholesky_cheatah(benchmark::State& state) {

197 const long long n = state.range(0);

198 const nd::NDArray a = cheatah_sym(n);

199 for (auto _ : state) benchmark::DoNotOptimize(la::cholesky(a));

200}

201void BM_cholesky_eigen(benchmark::State& state) {

202 const long long n = state.range(0);

203 const Eigen::MatrixXd a = eigen_sym(n);

204 for (auto _ : state) { Eigen::MatrixXd L = a.llt().matrixL(); benchmark::DoNotOptimize(L.data()); }

205}

206BENCHMARK(BM_cholesky_cheatah)->Arg(64);

207BENCHMARK(BM_cholesky_eigen)->Arg(64);

208

209// ---- qr -----------------------------------------------------------------

210void BM_qr_cheatah(benchmark::State& state) {

211 const long long n = state.range(0);

212 const nd::NDArray a = cheatah_matrix(n);

213 for (auto _ : state) benchmark::DoNotOptimize(la::qr(a).r);

214}

215void BM_qr_eigen(benchmark::State& state) {

216 const long long n = state.range(0);

217 const Eigen::MatrixXd a = eigen_matrix(n);

218 for (auto _ : state) {

219 Eigen::HouseholderQR<Eigen::MatrixXd> qr(a);

220 Eigen::MatrixXd r = qr.matrixQR().triangularView<Eigen::Upper>();

221 benchmark::DoNotOptimize(r.data());

222 }

223}

224BENCHMARK(BM_qr_cheatah)->Arg(32)->Arg(64);

225BENCHMARK(BM_qr_eigen)->Arg(32)->Arg(64);

226

227// ---- eigh (eigenvalues + eigenvectors) ----------------------------------

228void BM_eigh_cheatah(benchmark::State& state) {

229 const long long n = state.range(0);

230 const nd::NDArray a = cheatah_sym(n);

231 for (auto _ : state) benchmark::DoNotOptimize(la::eigh(a).values);

232}

233void BM_eigh_eigen(benchmark::State& state) {

234 const long long n = state.range(0);

235 const Eigen::MatrixXd a = eigen_sym(n);

236 for (auto _ : state) {

237 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(a); // values + vectors

238 benchmark::DoNotOptimize(es.eigenvalues().data());

239 }

240}

241BENCHMARK(BM_eigh_cheatah)->Arg(32)->Arg(64);

242BENCHMARK(BM_eigh_eigen)->Arg(32)->Arg(64);

243

244// ---- outer product ------------------------------------------------------

245void BM_outer_cheatah(benchmark::State& state) {

246 const long long n = state.range(0);

247 const nd::NDArray u = nd::full({n}, 1.5), v = nd::full({n}, 2.0);

248 for (auto _ : state) benchmark::DoNotOptimize(la::outer(u, v));

249}

250void BM_outer_eigen(benchmark::State& state) {

251 const long long n = state.range(0);

252 const Eigen::VectorXd u = Eigen::VectorXd::Constant(n, 1.5), v = Eigen::VectorXd::Constant(n, 2.0);

253 for (auto _ : state) { Eigen::MatrixXd c = u * v.transpose(); benchmark::DoNotOptimize(c.data()); }

254}

255BENCHMARK(BM_outer_cheatah)->Arg(64)->Arg(256);

256BENCHMARK(BM_outer_eigen)->Arg(64)->Arg(256);

257

258// ---- trace + Frobenius norm ---------------------------------------------

259void BM_trace_cheatah(benchmark::State& state) {

260 const nd::NDArray a = cheatah_matrix(256);

261 for (auto _ : state) benchmark::DoNotOptimize(la::trace(a));

262}

263void BM_trace_eigen(benchmark::State& state) {

264 const Eigen::MatrixXd a = eigen_matrix(256);

265 for (auto _ : state) benchmark::DoNotOptimize(a.trace());

266}

267BENCHMARK(BM_trace_cheatah);

268BENCHMARK(BM_trace_eigen);

269

270void BM_norm_cheatah(benchmark::State& state) {

271 const long long n = state.range(0);

272 const nd::NDArray a = cheatah_matrix(n);

273 for (auto _ : state) benchmark::DoNotOptimize(la::norm(a));

274}

275void BM_norm_eigen(benchmark::State& state) {

276 const long long n = state.range(0);

277 const Eigen::MatrixXd a = eigen_matrix(n);

278 for (auto _ : state) benchmark::DoNotOptimize(a.norm());

279}

280BENCHMARK(BM_norm_cheatah)->Arg(32)->Arg(256);

281BENCHMARK(BM_norm_eigen)->Arg(32)->Arg(256);

282

283} // namespace