  Most numerical methods for conic problems use the homogenous primal-dual
embedding, which yields a primal-dual solution or a certificate establishing
primal or dual infeasibility. Following Patrinos (and others, 2018), we express
the embedding as the problem of finding a zero of a mapping containing a
skew-symmetric linear function and projections onto cones and their duals. We
focus on the special case when this mapping is regular, i.e., differentiable
with nonsingular derivative matrix, at a solution point. While this is not
always the case, it is a very common occurrence in practice. We propose a
simple method that uses LSQR, a variant of conjugate gradients for least
squares problems, and the derivative of the residual mapping to refine an
approximate solution, i.e., to increase its accuracy. LSQR is a matrix-free
method, i.e., requires only the evaluation of the derivative mapping and its
adjoint, and so avoids forming or storing large matrices, which makes it
efficient even for cone problems in which the data matrices are given and
dense, and also allows the method to extend to cone programs in which the data
are given as abstract linear operators. Numerical examples show that the method
almost always improves an approximate solution of a conic program, and often
dramatically, at a computational cost that is typically small compared to the
cost of obtaining the original approximate solution. For completeness we
describe methods for computing the derivative of the projection onto the cones
commonly used in practice: nonnegative, second-order, semidefinite, and
exponential cones. The paper is accompanied by an open source implementation.
