When is the Rayleigh quotient (spherically) convex?
I provide a partial answer below, but I would be interested to see what people think.
Claim: If $A$ is not a multiple of the identity, then $f$ is not convex.
Let $\lambda_i$, $e_i$ denote the (real) eigenvalue-eigenvector pairs for $A$, for $i = 1, 2, \dots, n$.
Suppose that $\lambda_i < \lambda_j$ for some distinct $i, j \in \{1, 2, \dots, n\}$. Let us consider the geodesics with images lying in $E_{ij} := \mathrm{span}\{e_i, e_j\} \cap S^{n-1}$.
For $x \in E_{ij}$, write $x = c_i e_i + c_j e_j$, for $c_i^2 + c_j^2 = 1$, and note that $f|_{E_{ij}}(x) = \lambda_i c_i^2 + \lambda_j c_j^2$. We may reparameterize this as $f|_{E_{ij}}(\theta) = \lambda_i \cos^2(\theta) + \lambda_j \sin^2(\theta)$.
Geodesics $\gamma$ with images in $E_{ij}$ correspond to connected intervals $I \subset [0, 2\pi)$ (of length less than $\pi$). Hence such maps $f \circ \gamma$ are restrictions of $f_{E_{ij}}$ to intervals.
- There exists a connected interval $I$ of length less than $\pi$ on which $\cos^2(\theta)$ is neither convex nor concave.
- Hence, if $\lambda_i < \lambda_j$, $f \circ \gamma$ is not convex for some $\gamma$ whose image lies in $E_{ij}$. Thus $f$ is not convex.
- If no such pair $i, j$ exists, then $A = cI$ for some $c \in \mathbf{R}$.
The notion of convexity you are using is unnatural for functions defined on compact (connected) Riemannian manifolds: With this definition every convex function is constant. Indeed, suppose that $f: M\to {\mathbb R}$ is a convex function on a compact connected Riemannian manifold. Then $f$ is necessarily convex and, hence, attains its maximum at some point $p$. But for every $q\in M$ there is a geodesic $c$ in $M$ containing $q$ and $p$ as its interior point. Thus, the convex function $f|c$ is nonconstant. But a nonconstant convex function on an interval cannot attain its maximum at an interior point. Hence, $f(p)=f(q)$ and, therefore, $f$ is constant. (A small modification of this proof works even if you assume that $f$ is convex only on all distance-minimizing geodesics.)
Applied in the setting of your question, you have the bilinear form $\langle x, y\rangle=x^TAy$ on ${\mathbb R}^n$. Either the corresponding quadratic form $q(x)=f(x)$ is identically zero (which implies that $A=0$), or after multiplying the matrix $A$ by a scalar, we get $q(x)=1$ for all $x\in S^{n-1}$. But a bilinear form is uniquely determined by its quadratic form: $$ \langle x, y\rangle= \frac{1}{2}(q(x+y) - q(x) -q(y)). $$ Hence, $\langle x, y\rangle$ equals the standard dot product. In other words, $A=I$.