def conjugate_gradient_method(matrix_vector_operator, loss_gradient, nsteps, rdotr_tol=<NUM_LIT>):
x = torch.zeros_like(loss_gradient)<EOL>r = loss_gradient.clone()<EOL>p = loss_gradient.clone()<EOL>rdotr = torch.dot(r, r)<EOL>for i in range(nsteps):<EOL><INDENT>Avp = matrix_vector_operator(p)<EOL>alpha = rdotr / torch.dot(p, Avp)<EOL>x += alpha * p<EOL>r -= alpha * Avp<EOL>new_rdotr = torch.dot(r, r)<EOL>betta = new_rdotr / rdotr<EOL>p = r + betta * p<EOL>rdotr = new_rdotr<EOL>if rdotr < rdotr_tol:<EOL><INDENT>break<EOL><DEDENT><DEDENT>return x<EOL>
Conjugate gradient algorithm