Weird runtime behavior of scipy.linalg.solve_triangular and trtrs

I want to understand the time complexity of scipy.linalg.solve_triangular, which calls trtrs from LAPACK under the hood, so I wrote the following benchmarking script:

from time import perf_counter

import numpy as np
import scipy

def benchmark(p):
L = np.ones((p, p))
b = np.zeros((p))

start_time = perf_counter()
for _ in range(10):
scipy.linalg.solve_triangular(L, b,
lower=True,
unit_diagonal=True,
overwrite_b=True,
check_finite=False)
end_time = perf_counter()

return end_time - start_time

if __name__ == "__main__":
elapsed = benchmark(1<<8)
print(elapsed)
elapsed = benchmark(1<<10)
print(elapsed)
elapsed = benchmark(1<<12)
print(elapsed)
elapsed = benchmark(1<<14)
print(elapsed)
elapsed = benchmark(1<<16)
print(elapsed)


The output is

0.12597245816141367
0.03438423713669181
0.05818071821704507
0.1374253649264574
1.6546398717910051


This is surprising: why does solving a larger problem take less time, as in 1<<10 and 1<<12 vs 1<<8?

The more important question is: what's the time complexity of scipy.linalg.solve_triangular for $$L \in \mathbb{R}^{p \times p}$$ and $$b \in \mathbb{R}^{p}$$? I would expect it to be $$O(p^2)$$ with forward/backward substitution, but rumors are that triangular solve may take $$O(p^3)$$ in some implementations.

• It's probably just some spurious startup / compilation time for the first call. What happens if you start from 1<<6 instead? Or swap the order of the calls? Commented Feb 27 at 7:15
• @FedericoPoloni Thanks for the pointer! I think I have figured it out. See my answer below! Commented Feb 28 at 4:28

Thanks to Federico Poloni's comment, I have revised the script as below. It looks like scipy.linalg.solve_triangular has time complexity $$O(p^2)$$, because the execution time quadruples as the problem size doubles.

from time import perf_counter

import numpy as np
import scipy

def benchmark(p):
L = np.ones((p, p))
b = np.zeros((p))

# Warm up
scipy.linalg.solve_triangular(L, b,
lower=True,
unit_diagonal=True,
overwrite_b=True,
check_finite=False)

start_time = perf_counter()
for _ in range(64):
scipy.linalg.solve_triangular(L, b,
lower=True,
unit_diagonal=True,
overwrite_b=True,
check_finite=False)
end_time = perf_counter()

return end_time - start_time

if __name__ == "__main__":
last_elapsed = None
for i in range(16):
elapsed = benchmark(1<<i)
if last_elapsed is not None:
print(i, elapsed, elapsed/last_elapsed)
else:
print(i, elapsed)
last_elapsed = elapsed