0
$\begingroup$

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.

$\endgroup$
2
  • 1
    $\begingroup$ 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? $\endgroup$ Commented Feb 27 at 7:15
  • 1
    $\begingroup$ @FedericoPoloni Thanks for the pointer! I think I have figured it out. See my answer below! $\endgroup$
    – nalzok
    Commented Feb 28 at 4:28

1 Answer 1

4
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.