I'll begin with a general remark: first-order information (i.e., using only gradients, which encode slope) can only give you directional information: It can tell you that the function value decreases in the search direction, but not for how long. To decide how far to go along the search direction, you need extra information (gradient descent with constant step lengths can fail even for convex quadratic problems). For this, you basically have two choices:
- Use second-order information (which encodes curvature), for example by using Newton's method instead of gradient descent (for which you can always use step length $1$ sufficiently close to the minimizer).
- Trial and error (by which of course I mean using a proper line search such as Armijo).
If, as you write, you don't have access to second derivatives, and evaluating the obejctive function is very expensive, your only hope is to compromise: use enough approximate second-order information to get a good candidate step length such that a line search needs only $\mathcal{O}(1)$ evaluations (i.e., at most a (small) constant multiple of the effort you need to evaluate your gradient).
One possibility is to use Barzilai--Borwein step lengths (see, e.g., Fletcher: On the Barzilai-Borwein method. Optimization and control with applications, 235–256, Appl. Optim., 96, Springer, New York, 2005). The idea is to use a finite difference approximation of the curvature along the search direction to get an estimate of the step size. Specifically, choose $\alpha_0>0$ arbitrary, set $g^0:=\nabla f(x^0)$ and then for $k=0,...$:
- Set $s^k = -\alpha_k^{-1} g^k$ and $x^{k+1}=x^k+s^k$
- Evaluate $g^{k+1}=\nabla f(x^{k+1})$ and set $y^k = g^{k+1}-g^{k}$
- Set $\alpha_{k+1} = \frac{(y^k)^Ty^k}{(y^k)^Ts^k}$
This choice can be shown to converge (in practice very quickly) for quadratic functions, but the convergence is not monotone (i.e., the function value $f(x^{k+1})$ can be larger than $f(x^k)$, but only once in a while; see the plot on page 10 in Fletcher's paper). For non-quadratic functions, you need to combine this with a line search, which needs to be modified to deal with the non-monotonicity. One possibility is choosing $\sigma_k \in (0,\alpha_k^{-1})$ (e.g., by backtracking) such that
$$ f(x^k - \sigma_k g^k) \leq \max_{\max(k-M,1)\leq j\leq k} f(x^j) - \gamma \sigma_k (g^k)^Tg^k,$$
where $\gamma\in(0,1)$ is the typical Armijo parameter and $M$ controls the degree of monotonicity (e.g., $M=10$). There's also a variant that uses gradient values instead of function values, but in your case the gradient is even more expensive to evaluate than the function, so that doesn't make sense here. (Note: You can of course try to blindly accept the BB step lengths and trust your luck, but if you need any sort of robustness -- as you wrote in your comments -- that would be a really bad idea.)
An alternative (and, in my opinion, much better) approach would be to use this finite difference approximation already in the computation of the search direction; this is called a quasi-Newton method. The idea is to incrementally build an approximation of the Hessian $\nabla^2 f(x^k)$ by using differences of gradients. For example, you could take $H_0=\mathrm{Id}$ (the identity matrix) and for $k=0,\dots$ solve
$$H_{k}s^{k} = -g^{k},\label{cc1}\tag{1}$$
and set
$$H_{k+1} = H_k + \frac{(y^k-H_ks^k)^T(s^k)^T}{(s^k)^Ts^k}$$
with $y^k$ as above and $x^{k+1} = x^k +s^k$. (This is called Broyden update and is rarely used in practice; a better but slightly more complicated update is the BFGS update, for which -- and more information -- I refer to Nocedal and Wright's book Numerical Optimization.) The downside is that a) this would require solving a linear system in each step (but only of the size of the unknown which in your case is an initial condition, hence the effort should be dominated by solving PDEs to get the gradient; also, there exist update rules for approximations of the inverse Hessian, which only require computing a single matrix-vector product) and b) you still need a line search to guarantee convergence...
Luckily, in this context there exists an alternative approach that makes use of every function evaluation. The idea is that for $H_k$ symmetric and positive definite (which is guaranteed for the BFGS update), solving \eqref{cc1} is equivalent to minimizing the quadratic model
$$q_k(s) = \frac12 s^T H_k s + s^T g^k.$$
In a trust region method, you would do so with the additional constraint that $\|s\| \leq \Delta_k$, where $\Delta_k$ is an appropriately chosen trust region radius (which plays the role of the step length $\sigma_k$). The key idea is now to choose this radius adaptively, based on the computed step. Specifically, you look at the ratio
$$ \rho_k := \frac{f(x^k)-f(x^k+s^k)}{f(x^k)-q_k(s^k)}$$
of the actual and predicted reduction in function value. If $\rho_k$ is very small, your model was bad, and you discard $s^k$ and try again with $\Delta_{k+1}<\Delta_k$. If $\rho_k$ is close to $1$, your model is good, and you set $x^{k+1}=x^k+s^k$ and increase $\Delta_{k+1}>\Delta_k$. Otherwise you just set $x^{k+1}=x^k+s^k$ and leave $\Delta_k$ alone. To compute the actual minimizer $s^k$ of $\min_{\|s\|\leq \Delta_k} q_k(s)$, there exist several strategies to avoid having to solve the full constrained optimization problem; my favorite is Steihaug's truncated CG method. For more details, I again refer to Nocedal and Wright.
