#' Solve the matrix equation
#'
#'  X A X + X = B
#'
#' for symmetric, positive definite `A`, `B`. This is a special case of an
#' algebraic Riccati equation.
#'
#' @examples
#' A <- crossprod(rmatrix(5, 5))
#' B <- crossprod(rmatrix(5, 5))
#'
#' X <- riccati(A, B)
#'
#' all.equal(X %*% A %*% X + X, B)
#'
#' @export
riccati <- function(A, B, max.iter = 20L) {
    p <- nrow(A)

    V <- diag(0.5 - (p < seq_len(2 * p)))
    V[seq_len(p), seq_len(p) + p] <- A
    V[seq_len(p) + p, seq_len(p)] <- B

    for (iter in seq_len(max.iter)) {
        V <- 0.5 * (V + pinv(V))
    }

    G <- V + diag(2 * p)
    X <- pinv(G[, 1:p + p]) %*% -G[, 1:p]

    structure(X, resid = norm(X %*% A %*% X + X - B, "F"))
}