I have come across the following proof which proves the cyclic property of trace usually in matrix notation, however in this case, in bra-ket notation.
Imagine two operators $\hat{A}$ and $\hat{B}$, along with a Hermitian operator $\hat{F}$.
$\sum_i = \left<f_i\left|\hat{A}\hat{B}\right|f_i\right>$
$\sum_{i,j} = \left<f_i\left|\hat{A} \hat{n} \hat{B}\right|f_i\right>$, where the identity operator $\hat{n} = \sum_j \left|f_j\left>\right<f_j\right|$,
$\sum_{i,j} = \left<f_i\left|\hat{A} \right|f_j\right>\left<f_j\left| \hat{B}\right|f_i\right>$
$\sum_{i,j} =\left<f_j\left| \hat{B}\right|f_i\right>\left<f_i\left|\hat{A} \right|f_j\right>$.
$\vdots$
I am unsure however, which property was used such that $\left<f_j\left| \hat{B}\right|f_i\right>$ and $\left<f_i\left|\hat{A} \right|f_j\right>$ were interchanged.
This is because these are "scalar products" that are $\mathbb{C}$-valued, and multiplication is commutative in $\mathbb{C}$.
In case $A$, $B$, $F$ are matrices, $\big\langle f_i\big|A\big|f_j\big\rangle$ is the same as $\overline{f_i}^\text{T}Af_j\in\mathbb{C}$.
Note that in general, the operators do not have to be matrices, and the underlying Hilbert space can be an infinite-dimensional space of functions, such as $L^2[0,1]$ or $\ell^2$.