Asymptotic Complexity of b.transpose()*A^(-1)*b

Question

double sumtrv1(const Eigen::MatrixXd &A, const Eigen::VectorXd &b) {
    const int n = A.cols();
    assert((A.rows() == n) && (b.size() == n));

    return b.transpose() *
           A.triangularView<Eigen::Upper()>.solve(
               Eigen::MatrixXd::Identity(n,n)) *
           b;
}

so apparently, this code has asymptotic complexity of O(n^3) with the reasoning:

The asymptotic complexity is O(n^3), because the code is solving $n$ linear systems of equations with an (n x n) upper triangular system matrix. This amounts to $n$ backward substitutions, each of which costs O(n^2) operations.

My thought would have been: -Transposing might cost n operations. Although we probably can neglect it since we just switch from row to col major or vice versa. (No idea what Eigen really does here) We then get the upper triangular part of A. That's also "free". We then get the inverse of this upper triangular part of A. This costs us $n^2$ We then basically calculate b.transpose()*A^(-1)b which is again nn

So overall, we get O(n*n^2*n)=O(n4)

What's wrong with my reasoning?

Answer 1

Yes, the complexity of this code is O(n^3) because computing the inverse of the triangular matrix A is O(n^3): you have n backward substitutions (one for each of the n columns of the identity matrix), and each backward substitution cost O(n^2). Then you have in addition one matrix vector product costing O(n^2), and one inner product costing O(n), but those two operations are negligible compared to O(n^3) and the overall complexity is O(n^3).

That being said, you can rewrite your expression to get a O(n^2) complexity:

b' * triu(A)^-1 * I * b = b' * triu(A)^-1 * b

In Eigen language, you get:

return b.transpose() * A.triangularView<Eigen::Upper()>.solve(b);

You now have one backward substitution O(n^2) plus one inner product O(n), yielding a O(n^2) asymptotic complexity. No need to say this new version will be considerably faster!