Arguably LU &#x2F; LDLT factorizations before solving the direct system are the more elegant &#x2F; beautiful &#x2F; simpler answers if you take the algorithm at face value. However, large-dimension matrix inversion has problems outside of just beauty. Part of the problem is that we don&#x27;t have fast, perfect representations of real or complex numbers. At best we have double precision floating points, which means that when you take the adjoint over the determinant, you&#x27;re looking at a lot of error.<p>That&#x27;s not to say that the implementations of these have to be ugly or that our current BLAS &#x2F; LAPACK &#x2F; whatever can&#x27;t be improved. There&#x27;s a lot of room for improvement in this space, but the fact that they&#x27;re 1. fast and 2. well tested and robust means that improving these sorts of interfaces come with baggage. People complain about the wheel being reinvented every time somebody invents a new Javascript framework, but could you imagine if this same problem applied to BLAS and other high performance math libraries? (Arguably, at least then we&#x27;d have more choices in math space)