Disclaimer: This follows lecture notes by Prof. Reusken at RWTH Aachen University who cites W. Hackbusch. Iterative Solution of Large Sparse Systems of Equations as the main resource.
There are most likely better motivations for setting up the function $$f(x) := \boldsymbol x^T A \boldsymbol x + \boldsymbol b^T \boldsymbol x = \langle \boldsymbol x, A \boldsymbol x\rangle - \langle \boldsymbol x, \boldsymbol b \rangle \tag1$$
but I think you can argue a bit like this:
Okay, let's say we want to solve the linear system $$A\boldsymbol x=\boldsymbol b \tag2$$ with symmetric positive matrix (s.p.d) $A$. From very basic multidimensional optimization, one knows that if the Hessian $\nabla^2 f(\boldsymbol x) $ of the objective function $f(\boldsymbol x)$ is positive definite $\forall \: \boldsymbol x$, the objective function is globally strictly convex and thus there is only one global minimizer.
Let's try to leverage the s.p.d. property of $A$: We aim for an objective function such that $$\nabla^2 f(\boldsymbol x) \overset{!}{=} A. \tag3$$
Again from optimization theory for convex functions, we know that we are at the optimum $\boldsymbol x^\star$ if $\nabla f(\boldsymbol x^\star) = \boldsymbol 0 $. The optimal $\boldsymbol x^\star$ for our problem would be the one that solves the linear system: $$A\boldsymbol x^\star =\boldsymbol b.\tag4$$
We see that this can be enforced by "integrating/solving the PDE" $\nabla^2 f(\boldsymbol x) = A$ and picking the integration constant in a clever way, namely as $-\boldsymbol b$:
$$\nabla f(x) = A \boldsymbol x - \boldsymbol b. \tag5$$
Again, by clever integration / PDE solve, you can construct the objective function $$f(\boldsymbol x) = \frac12 \boldsymbol x^T A \boldsymbol x - \boldsymbol x^T \boldsymbol b. \tag6$$
Now solving the linear system is equivalent to minimizing $f(\boldsymbol x)$. Again, this is just a motivation. Other ways of arriving at $f(\boldsymbol x)$ could stem from experience or constructing the sort of simplest convex function involving $A, \boldsymbol b$.
Now continue with optimization 101-stuff: The first method you learn in numerical optimization is steepest descent, which can be slightly more generalized written as
$$ \boldsymbol x^{k+1} = \boldsymbol x^k + \alpha \boldsymbol p^{k}. \label{7} \tag7$$
In the case of steepest descent, $\boldsymbol p^k = -\nabla f (\boldsymbol x^k)$. Let us leave $\boldsymbol p^k \neq \boldsymbol 0 $ general for the moment and let's turn our attention towards the step-length $\alpha$:
$$f\big(\boldsymbol x^{k+1} \big) = f\big(\boldsymbol x^k + \alpha \boldsymbol p^{k}\big) = f\big(\boldsymbol x^k \big) + \alpha \big\langle \boldsymbol p^k , \underbrace{A \boldsymbol x^k -\boldsymbol b}_{=:\boldsymbol r^k} \big\rangle + \frac12 \alpha^2 \big\langle \boldsymbol p^k , A \boldsymbol p^k \big\rangle =: g(\alpha) \label{8}\tag8$$
For fixed $\boldsymbol x^k, \boldsymbol p^k$ this is a scalar, convex, univariate function in $\alpha$ with minimum at $\alpha^\star$ such that $g'(\alpha^\star) = 0$ which is clearly the case when $$\alpha^\star \big(\boldsymbol p^k, \boldsymbol r^k \big) =- \frac{ \big \langle \boldsymbol p^k, \boldsymbol r^k \big \rangle }{\big \langle \boldsymbol p^k,A \boldsymbol p^k \big \rangle}. \label{9} \tag9$$
Again, valid for any choice of search-direction $ \boldsymbol p^k$!
It is instructive to reconsider the steepest descent method. At a point $\boldsymbol x ^k \neq \boldsymbol x^\star $ the gradient $\nabla f\big(\boldsymbol x^k \big)$ is given by
$$\nabla f\big(\boldsymbol x^k \big) = A \boldsymbol x^k - \boldsymbol b = \boldsymbol r^k \label{10} \tag{10}$$
where we call $\boldsymbol r^k$ the residual of the k'th iterate, since it quantifies the defect/error of the current iterate and the desired outcome. Besides the intuitive reasoning of picking $\boldsymbol p^k$ as the negative gradient, there is another one available. Recall \eqref{8} and note that $$g'(\alpha) = \frac{\partial}{\partial \alpha} f\big(\underbrace{\boldsymbol x^k + \alpha \boldsymbol p^{k}}_{=:\boldsymbol z}\big) = \nabla f(\boldsymbol y) \frac{\partial}{\partial \alpha} \boldsymbol z = \nabla f(\boldsymbol z) \cdot \boldsymbol p^k \tag{11}$$ which is for $\alpha = 0$ the precise definition of the directional derivative of $f\big(\boldsymbol x^k\big)$ in direction $\boldsymbol p^k$! The minimum among all search directions based on knowledge of only $\boldsymbol x^k$ is thus the one for which $\nabla f(\boldsymbol y) \cdot \boldsymbol p^k = g'(0) $ is minimized. Here we actually maximize the absolute value $ |g'(0) |$ since the nominator of $\alpha$ \eqref{9} takes care of selecting the correct sign.
$$ \max_{\boldsymbol p^k} \Big | \big\langle \boldsymbol p^k , \underbrace{A \boldsymbol x^k -b }_{= \boldsymbol r^k}\big \rangle \Big | \tag{12}$$
and thus $$\boldsymbol p^k = \boldsymbol r^k = A \boldsymbol x^k - \boldsymbol b = \nabla f\big(\boldsymbol x^k \big). \label{13} \tag{13}$$
Sometimes you require that $\boldsymbol p^k$ is normalized, i.e., $ \big \Vert \boldsymbol p^k \big \Vert_2 = 1$, but again, all of this is basically covered by $\alpha$.
Again, the important part is that the search direction is optimal w.r.t. the current iterate $\boldsymbol x^k$.
It is instructive to consider also the relation of two successive search directions $\boldsymbol p^k = \boldsymbol r^k, \boldsymbol p^{k+1} = \boldsymbol r^{k+1}$.
First, observe that from \eqref{7} you can deduce that $$\boldsymbol r^{k+1} = A \boldsymbol x^{k+1} - \boldsymbol b \overset{\eqref{7}}{=} A \big(\boldsymbol x^{k} + \alpha^\star \boldsymbol p^{k} \big) - \boldsymbol b = \boldsymbol r^{k} + \alpha^\star A \boldsymbol p^{k} \overset{\eqref{13}}{=} \boldsymbol r^{k} + \alpha^\star A \boldsymbol r^{k} \label{14}\tag{14}$$
Then examine \begin{align} \big \langle \underbrace{\boldsymbol p^k}_{ = \boldsymbol r^k}, \underbrace{\boldsymbol p^{k+1}}_{ = \boldsymbol r^{k+1}} \big \rangle &= \big \langle \boldsymbol r^k, \boldsymbol r^{k} + \alpha^\star A \boldsymbol r^{k} \big \rangle \tag{15a}\\
&= \big \langle \boldsymbol r^k, \boldsymbol r^k \big \rangle + \alpha^\star \big \langle \boldsymbol r^k, A \boldsymbol r^{k} \big \rangle \tag{15b}\\
&\overset{\eqref{9}, \eqref{13}}{=} \big \langle \boldsymbol r^k, \boldsymbol r^k \big \rangle - \underbrace{\frac{ \big \langle \boldsymbol r^k, \boldsymbol r^k \big \rangle }{\big \langle \boldsymbol r^k,A \boldsymbol r^k \big \rangle}}_{\alpha^\star} \big \langle \boldsymbol r^k, A\boldsymbol r^k \big \rangle = 0 \label{15c}\tag{15c}\end{align}
So two successive search directions are orthogonal w.r.t. the standard/Euclidian scalar product! Orthogonality is also something that shows up for Conjugate Gradients, but then w.r.t. a different scalar product.
To be able to finally proceed to CG, we first need a definition:
For $V$ being a subspace of $\mathbb R^n$, i.e., $V \subset \mathbb R^n$ we call $\boldsymbol y \in \mathbb R^n$ optimal for $V$ if $$f(\boldsymbol y) = \min_{\boldsymbol u \in V} f(\boldsymbol y + \boldsymbol u). \label{16}\tag{16}$$
Now let $\big \{ \boldsymbol d_i \big \}_{i = 1, \dots s} $ be a basis of $V$. Furthermore, for given $\boldsymbol y$ and $\boldsymbol c \in \mathbb R^s $ define $$h(\boldsymbol c) := f\bigg(\boldsymbol y + \sum_{i=1}^s c_i \boldsymbol d_i \bigg) \tag{17}$$
Then $\boldsymbol y$ is optimal for $V$ iff $\nabla_{\boldsymbol c} h(\boldsymbol 0) = \boldsymbol 0 $. This follows from the definition \eqref{16}. Carrying out the the differentiation, one obtains by employing the Chain Rule the "condition for optimality for subspaces"
$$ \nabla_{\boldsymbol c} h(\boldsymbol 0) = \big \langle f(\boldsymbol y), \boldsymbol d_i \big \rangle = \big \langle A \boldsymbol y - \boldsymbol b, \boldsymbol d_i \big \rangle \overset{\eqref{10}}{=} \big \langle \boldsymbol r_{\boldsymbol y}, \boldsymbol d_i \big \rangle\overset{!}{=} 0 \label{18}\tag{18} $$
From \eqref{15c} and \eqref{14} it is clear that the steepest descent iterate $\boldsymbol x^{k+1}$ is optimal to the subspace $\text{span} \big \{ \boldsymbol p^k \big \} = \text{span} \big \{ \boldsymbol r^k \big \}$.
This can serve as the motivation to find now iterates that are not only optimal w.r.t. the subspace spanned by the very previous search direction, but all previous search directions. By considering some small (2D) examples, one can see that the one can indeed find such iterates, which then lead to the (later generalized) property, namely that the CG method reaches the analytical solution in exactly $n$ steps for a $n$-dimensional problem.
The following derivation is induction-like: Use the steepest descent method to initialize, i.e., $\boldsymbol p^0 = \boldsymbol r^0$ and $\boldsymbol x^1$ computed by \eqref{7}. Then, $\boldsymbol x^1$ is optimal for $\text{span} \big \{ \boldsymbol p^0 \big \}$. By standard induction, assume now that for a certain iteration $k$, $1\leq k < n$, the current iterate $\boldsymbol x^k$ is optimal for $V_k := \text{span} \big \{ \boldsymbol p^0, \dots \boldsymbol p^{k-1} \big \}$. Also assume that we are not yet at the solution, thus $\boldsymbol r^k \neq \boldsymbol 0$. Now we want to pick $\boldsymbol p^k, \boldsymbol x^{k+1}$ such that $\boldsymbol x^{k+1}$ is optimal for $V_{k+1}$. It turns out that it is an excellent idea to choose $\boldsymbol p^k$ such that it it is orthogonal w.r.t. the $A$-scalar product $\langle \boldsymbol u, \boldsymbol v \rangle_A := \langle \boldsymbol u, A \boldsymbol v \rangle $, i.e.,
$$ \boldsymbol p^k \perp_A V_k \Leftrightarrow \boldsymbol p^k \in V_k^{\perp_A}. \tag{19}$$
This alone does not yet define a unique search direction $\boldsymbol p^k$. A unique one can be found by picking $\boldsymbol p^k \in V_k^{\perp_A}$ with maximum descent in $A$ norm $\Vert \boldsymbol u \Vert_A := \sqrt{\langle \boldsymbol u, \boldsymbol u \rangle_A}$ i.e.,
$$ \boldsymbol p^k = \text{argmin}_{\boldsymbol p \in V_k^{\perp_A}} \Vert \boldsymbol p - \boldsymbol r^k \Vert_A \tag{20}$$
This choice corresponds to $ \boldsymbol p^k$ being the $A$-orthogonal projection of $\boldsymbol r^k$ onto $V_k^{\perp_A}$.
Thus, $ \boldsymbol p^k$ can be determined by the Gram-Schmidt Orthogonalization/Projection process:
$$\boldsymbol p^k = \boldsymbol r^k - \sum_{i=0}^{k-1} \frac{ \big \langle \boldsymbol p^i, \boldsymbol r^k \big \rangle_A}{\big \langle \boldsymbol p^i, \boldsymbol p^i \big \rangle_A } \boldsymbol p^i. \label{21} \tag{21}$$
It remains to show that $\boldsymbol x^{k+1}$ is indeed optimal for $V_{k+1}$ for this choice of $\boldsymbol p^k$. By assumption $\boldsymbol x^k$ is optimal for $V_k$, thus $\langle \boldsymbol p^i, \boldsymbol r^k \rangle = 0$ for $i = 0, \dots k-1$. This implies $\boldsymbol r^k \perp V_k $. Since $\boldsymbol r^k \neq \boldsymbol 0 \Rightarrow \boldsymbol r^k \notin V_k $ and thus by \eqref{21} also $\boldsymbol p^k \notin V_k $. In other words, $\boldsymbol p^k$ is linearly independent of the previous $\boldsymbol p^i$ and $V_{k+1}$ has indeed dimension $k+1$. From equations \eqref{7}, \eqref{9} it follows that $$ \big \langle \boldsymbol p^k, - \boldsymbol r^{k+1} \big \rangle = \big \langle \boldsymbol p^k, - \boldsymbol r^{k} \big \rangle - \alpha^\star \big \langle \boldsymbol p^k, A \boldsymbol p^{k} \big \rangle = 0 \tag{22}.$$
For $1 \leq i < k$ we have from \eqref{7}
$$ \big \langle \boldsymbol p^i, - \boldsymbol r^{k+1} \big \rangle \overset{\eqref{10}}{=} \big \langle \boldsymbol p^i, \boldsymbol b - A \boldsymbol x^{k+1} \big \rangle \overset{\eqref{7}}{=} \big \langle \boldsymbol p^i, \boldsymbol b - A \boldsymbol x^k \big \rangle - \alpha^\star \langle \boldsymbol p^i, A \boldsymbol p^k \big \rangle \tag{23}.$$ where the first scalar product is zero due to the optimality of $\boldsymbol x^k $ w.r.t subspace $V_k$ and the second scalar product is zero due to the definition of $\alpha^\star$, \eqref{9}.
Now some I will show some important practical properties. Especially the Gram-Schmidt projection process \eqref{21} can be significantly simplified - this leads to then to the familiar $\beta^k$ of the CG Algortihm and a simpler update formula for $\boldsymbol p^k$.
We need for this the fact that $V_k := \text{span} \{ \boldsymbol p^0, \dots, \boldsymbol p^{k-1} \} = \text{span} \{ \boldsymbol r^0, \dots, \boldsymbol r^{k-1} \}. =: R_k$
This can be proven by induction. The hypothesis for $k = 0$ trivially holds since the first step done is the steepest descent one with $V_0 = \boldsymbol p^0 = \boldsymbol r^0 = R_0$. For $V_{k+1} = \text{span}\{V_k, \boldsymbol p^k \}$ we have from hypothesis that $V_k = R_k$. Thus, we have to show that $\boldsymbol p^k \in R_{k+1}$. \eqref{21} immediately shows that $\boldsymbol p^k \in \text{ span} \{\boldsymbol p^0, \dots, \boldsymbol p^{k-1}, \boldsymbol r^k \}$ which is by induction hypothesis identical to $\boldsymbol p^k \in \text{ span} \{R_k, \boldsymbol r^k \} = R_{k+1}.$ This shows the claim.
To finally have the familiar update formula for $\boldsymbol p^k$ we recall that $\boldsymbol x^k$ is optimal for the subspace $V_k$. By \eqref{18}, this is formalized as
$\big \langle \boldsymbol r^k, \boldsymbol p^j \big \rangle = 0 , j = 0, \dots, k-1$.
It was just shown that $V_k = R_k = \text{span}\{\boldsymbol r^0, \dots, \boldsymbol r^{k-1} \}$ and thus we have that $$\big \langle \boldsymbol r^k, \boldsymbol r^j \big \rangle = 0, j = 0, \dots, k -1 \tag{24}$$
Recall now the recursion for the residuals, \eqref{14}. From this formula it is obvious that $A \boldsymbol p ^j \in \text{span} \{\boldsymbol r^k,\boldsymbol r^{k +1}\}$. As a consequence,
$$ \big \langle \boldsymbol p^i, \boldsymbol r^k \big \rangle_A = \big \langle \boldsymbol Ap^i, \boldsymbol r^k \big \rangle = 0, i = 0, \dots, k-2 \tag{25}$$ and $$\boldsymbol p^k = \boldsymbol r^k - \frac{ \big \langle \boldsymbol p^{k-1}, \boldsymbol r^k \big \rangle_A}{\big \langle \boldsymbol p^{k-1}, \boldsymbol p^{k-1} \big \rangle_A } \boldsymbol p^{k-1}. \tag{26}$$
This finally gives the familiar algorithm:
$\boldsymbol x^0 $ is a given starting vector, compute initial residual $\boldsymbol r^0 = A \boldsymbol x^0 - \boldsymbol b = \boldsymbol p^0$, k = 0.
while ($\boldsymbol r^0 \neq \boldsymbol 0$ && k < k_max) {
k++; $\boldsymbol p^k = \boldsymbol r^k - \frac{ \big \langle \boldsymbol r^{k}, A\boldsymbol p^{k-1} \big \rangle}{\big \langle \boldsymbol p^{k-1}, A\boldsymbol p^{k-1} \big \rangle } \boldsymbol p^{k-1}$
$ \boldsymbol x^{k+1} = \boldsymbol x^k - \frac{ \big \langle \boldsymbol p^k, \boldsymbol r^k \big \rangle }{\big \langle \boldsymbol p^k,A \boldsymbol p^k \big \rangle} \boldsymbol p^{k} $
$\boldsymbol r^{k+1} =\boldsymbol r^{k} + \frac{ \big \langle \boldsymbol p^k, \boldsymbol r^k \big \rangle }{\big \langle \boldsymbol p^k,A \boldsymbol p^k \big \rangle} A \boldsymbol p^{k}$
}
