Many solutions to constrained optimization problems require finding the solution to a linear system Ax = b. It is commonly the case that solving a linear system is the most costly step in optimizing the objective function and consequently, much work has been done in speeding up various solvers by exploiting various properties or structures within the matrix. For example, one can use the LU or QR decomposition to solve general linear systems of the form Ax = b. However, if the matrix A satisfies certain properties, then often the computation can be sped up by using a different matrix factorization. For example, for positive definite matrices there are solvers using the Cholesky decomposition, or the for more general symmetric matrices there are solvers using the LDLT decomposition. Furthermore, more recent work has focused on optimizing solvers for Ax = b where A is a sparse matrix, resulting in sparse versions of the aforementioned solvers.