The following algorithm computes the Least squares solution || Ax - b|| subject to the equality constrain Bx = d. It's a classic algorithm that can be implemented only using a QR decomposition and a least squares solver. This implementation uses numpy and scipy. It makes use of the new linalg.solve_triangular function …