def conjugate_gradient_normal(op, x, rhs, niter=1, callback=None):
    """Optimized implementation of CG for the normal equation.

    This method solves the inverse problem (of the first kind) ::

        A(x) == rhs

    with a linear `Operator` ``A`` by looking at the normal equation ::

        A.adjoint(A(x)) == A.adjoint(rhs)

    It uses a minimum amount of memory copies by applying re-usable
    temporaries and in-place evaluation.

    The method is described (for linear systems) in a
    `Wikipedia article
    <https://en.wikipedia.org/wiki/Conjugate_gradient_method#\
Conjugate_gradient_on_the_normal_equations>`_.

    Parameters
    ----------
    op : `Operator`
        Operator in the inverse problem. If not linear, it must have
        an implementation of `Operator.derivative`, which
        in turn must implement `Operator.adjoint`, i.e.
        the call ``op.derivative(x).adjoint`` must be valid.
    x : ``op.domain`` element
        Element to which the result is written. Its initial value is
        used as starting point of the iteration, and its values are
        updated in each iteration step.
    rhs : ``op.range`` element
        Right-hand side of the equation defining the inverse problem
    niter : int
        Number of iterations.
    callback : callable, optional
        Object executing code per iteration, e.g. plotting each iterate.

    See Also
    --------
    conjugate_gradient : Optimized solver for symmetric matrices
    odl.solvers.smooth.nonlinear_cg.conjugate_gradient_nonlinear :
        Equivalent solver for the nonlinear case
    """
    # TODO: add a book reference
    # TODO: update doc

    if x not in op.domain:
        raise TypeError('`x` {!r} is not in the domain of `op` {!r}'
                        ''.format(x, op.domain))

    d = op(x)
    d.lincomb(1, rhs, -1, d)               # d = rhs - A x
    p = op.derivative(x).adjoint(d)
    s = p.copy()
    q = op.range.element()
    sqnorm_s_old = s.norm() ** 2  # Only recalculate norm after update

    for _ in range(niter):
        op(p, out=q)                       # q = A p
        sqnorm_q = q.norm() ** 2
        if sqnorm_q == 0.0:  # Return if residual is 0
            return

        a = sqnorm_s_old / sqnorm_q
        x.lincomb(1, x, a, p)               # x = x + a*p
        d.lincomb(1, d, -a, q)              # d = d - a*Ap
        op.derivative(p).adjoint(d, out=s)  # s = A^T d

        sqnorm_s_new = s.norm() ** 2
        b = sqnorm_s_new / sqnorm_s_old
        sqnorm_s_old = sqnorm_s_new

        p.lincomb(1, s, b, p)               # p = s + b * p

        if callback is not None:
            callback(x)
Optimized implementation of CG for the normal equation.

    This method solves the inverse problem (of the first kind) ::

        A(x) == rhs

    with a linear `Operator` ``A`` by looking at the normal equation ::

        A.adjoint(A(x)) == A.adjoint(rhs)

    It uses a minimum amount of memory copies by applying re-usable
    temporaries and in-place evaluation.

    The method is described (for linear systems) in a
    `Wikipedia article
    <https://en.wikipedia.org/wiki/Conjugate_gradient_method#\
Conjugate_gradient_on_the_normal_equations>`_.

    Parameters
    ----------
    op : `Operator`
        Operator in the inverse problem. If not linear, it must have
        an implementation of `Operator.derivative`, which
        in turn must implement `Operator.adjoint`, i.e.
        the call ``op.derivative(x).adjoint`` must be valid.
    x : ``op.domain`` element
        Element to which the result is written. Its initial value is
        used as starting point of the iteration, and its values are
        updated in each iteration step.
    rhs : ``op.range`` element
        Right-hand side of the equation defining the inverse problem
    niter : int
        Number of iterations.
    callback : callable, optional
        Object executing code per iteration, e.g. plotting each iterate.

    See Also
    --------
    conjugate_gradient : Optimized solver for symmetric matrices
    odl.solvers.smooth.nonlinear_cg.conjugate_gradient_nonlinear :
        Equivalent solver for the nonlinear case