How to understand the trace operation below?

I think that by $''$cyclic property of the trace$''$ we mean the following :

Consider $:k:$ finite dimensional matrices begin{align} mathrm A_1 & boldsymbol{=} n_1 times n_2 quad texttt{matrix} tag{1-01}label{1-01} mathrm A_2 & boldsymbol{=} n_2 times n_3 quad texttt{matrix} tag{1-02}label{1-02} mathrm A_3 & boldsymbol{=} n_3 times n_4 quad texttt{matrix} tag{1-03}label{1-03} cdots & boldsymbol{=} cdotcdot: times cdotcdot: quad texttt{matrix} tag{1-....}label{1-....} mathrm A_k & boldsymbol{=} n_k times n_1 quad texttt{matrix} tag{1-0k}label{1-0k} end{align} Then the product matrix $:mathrm A_1mathrm A_2mathrm A_3cdotsmathrm A_k:$ and all cyclic permutations of it are square matrices. Explicitly begin{align} mathrm A_1mathrm A_2mathrm A_3cdotsmathrm A_k & boldsymbol{=} n_1 times n_1 quad texttt{square matrix} tag{2-01}label{02-01} mathrm A_2mathrm A_3mathrm A_4cdotsmathrm A_1 & boldsymbol{=} n_2 times n_2 quad texttt{square matrix} tag{2-02}label{02-02} mathrm A_3mathrm A_4mathrm A_5cdotsmathrm A_2 & boldsymbol{=} n_3 times n_3 quad texttt{square matrix} tag{2-03}label{01-03} cdots & boldsymbol{=} cdotcdot: times cdotcdot: quad texttt{square matrix} tag{2-....}label{01-....} mathrm A_{rho}mathrm A_{rhoboldsymbol{+}1}mathrm A_{rhoboldsymbol{+}2}cdotsmathrm A_{rhoboldsymbol{-}1} & boldsymbol{=} n_rho times n_rho quad texttt{square matrix} tag{2-0k}label{01-0k} end{align}

Then begin{equation} texttt{Tr}left(mathrm A_1mathrm A_2mathrm A_3cdotsmathrm A_kright)boldsymbol{=}texttt{Tr}left(mathrm A_2mathrm A_3mathrm A_4cdotsmathrm A_1right)boldsymbol{=}texttt{Tr}left(mathrm A_3mathrm A_4mathrm A_5cdotsmathrm A_2right)boldsymbol{=}cdots tag{3}label{3} end{equation} Under this spirit if begin{align} mathrm A_1 & boldsymbol{equiv} langlepsi|boldsymbol{=} 1 times n boldsymbol{=} texttt{one row matrix} tag{4-01}label{04-01} mathrm A_2 & boldsymbol{equiv} A boldsymbol{=} n times n quad texttt{square matrix} tag{4-02}label{04-02} mathrm A_3 & boldsymbol{equiv} |psirangleboldsymbol{=} n times 1 boldsymbol{=} texttt{one column matrix} tag{4-03}label{04-03} end{align} then
begin{equation} texttt{Tr}left(mathrm A_1mathrm A_2mathrm A_3right)boldsymbol{=}texttt{Tr}left(mathrm A_2mathrm A_3mathrm A_1right)boldsymbol{=}texttt{Tr}left(mathrm A_3mathrm A_1mathrm A_2right) tag{5}label{5} end{equation} is translated to begin{equation} texttt{Tr}left(mathrm langlepsi|Amathrm |psirangleright)boldsymbol{=}texttt{Tr}left(Amathrm |psiranglemathrm langlepsi|right)boldsymbol{=}texttt{Tr}left(mathrm |psiranglemathrm langlepsi|Aright) tag{6}label{6} end{equation} But in any case is not permissible to write begin{equation} texttt{Tr}left(mathrm A_1mathrm A_2mathrm A_3right)boldsymbol{=}texttt{Tr}left(mathrm A_1right)texttt{Tr}left(mathrm A_2right)texttt{Tr}left(mathrm A_3right) tag{7}label{7} qquad textbf{(wrong !!!)} end{equation} as OP does in the question since on one hand the trace of non-square matrices has no sense and on the other hand this is wrong even in the case that all matrices are square.

Now, I think that a proof of the validity of the cyclic property eqref{3} for finite dimensional matrices $:rm A_rho, rho=1,2,cdots k:$ will help us to sketch a proof as an answer to the question. It's easy to realize that it's sufficient to prove eqref{3} for two matrices, say $:rm P={rm p_{ij}}:$ and $:rm Q={q_{kell}}:$ of dimensions $:ntimes m:$ and $:mtimes n:$ respectively begin{equation} texttt{Tr}left(mathrm Pmathrm Qright)boldsymbol{=}texttt{Tr}left(mathrm Qmathrm Pright) tag{8}label{8} end{equation} since for any number $:k:$ of matrices begin{equation} texttt{Tr}left(underbrace{mathrm A_1}_{rm P}underbrace{mathrm A_2mathrm A_3cdotsmathrm A_k}_{rm Q}right)boldsymbol{=}texttt{Tr}left(underbrace{mathrm A_2mathrm A_3cdotsmathrm A_k}_{rm Q}underbrace{mathrm A_1}_{rm P}right) tag{9}label{9} end{equation} The proof of eqref{8} runs as follows begin{equation} texttt{Tr}left(mathrm Pmathrm Qright)boldsymbol{=}sumlimits_{i=1}^{i=n}sumlimits_{j=1}^{j=m}mathrm p_{ij}mathrm q_{ji}boldsymbol{=}sumlimits_{j=1}^{j=m}sumlimits_{i=1}^{i=n}mathrm q_{ji}mathrm p_{ij}boldsymbol{=}texttt{Tr}left(mathrm Qmathrm Pright) tag{10}label{10} end{equation}

We sketch now a quick proof of equation eqref{6} in the special case of the 3-dimensional Hilbert space $:mathbb C^{3}$. Later on this proof will be generalized for an infinite-dimensional separable Hilbert space (with infinite countable complete orthonormal basis $:{phi_mu}$).

For the state vector $:|psirangle:$ we have begin{align} |psirangle & boldsymbol{=} z_1,phi_1boldsymbol{+}z_2,phi_2boldsymbol{+}z_3,phi_3boldsymbol{=}sumlimits_{mu=1}^{mu=3}z_mu,phi_mu,,quad z_muboldsymbol{=}langlephi_mu|psirangle in mathbb C tag{11a}label{11a} & texttt{so} quadquad |psirangle boldsymbol{=} sumlimits_{mu=1}^{mu=3}|phi_muranglelanglephi_mu|psirangle tag{11b}label{11b} langlepsi| & boldsymbol{=} overline{z_1},phi_1boldsymbol{+}overline{z_2},phi_2boldsymbol{+}overline{z_3},phi_3boldsymbol{=}sumlimits_{mu=1}^{mu=3}overline{z_mu},phi_mu,,quad overline{z_mu}boldsymbol{=}overline{langlephi_mu|psirangle}boldsymbol{=}langlepsi|phi_murangle in mathbb C tag{11c}label{11c} & texttt{so} quadquad langlepsi|boldsymbol{=} sumlimits_{mu=1}^{mu=3}langlephi_mu|langlepsi|phi_murangleboldsymbol{=}sumlimits_{mu=1}^{mu=3}langlephi_mu|overline{langlephi_mu|psirangle} tag{11d}label{11d} end{align} Formally we could consider them as one-column and one-row matrices respectively begin{equation} |psirangle boldsymbol{=} begin{bmatrix} :z_1 : vphantom{dfrac{a}{b}}: :z_2 : vphantom{dfrac{a}{b}}: :z_3 : vphantom{dfrac{a}{b}}: end{bmatrix},,qquad langlepsi| boldsymbol{=} begin{bmatrix} :overline{z_1} : & :overline{z_2} : & :overline{z_3} : vphantom{dfrac{a}{b}}: end{bmatrix} tag{12}label{12} end{equation} Consider that the operator $:A:$ is represented by the following $:3times 3:$ complex matrix begin{equation} Aboldsymbol{=} begin{bmatrix} a_{11} & a_{12} & a_{13} vphantom{dfrac{a}{b}}: a_{21} & a_{22} & a_{23} vphantom{dfrac{a}{b}}: a_{31} & a_{32} & a_{33} vphantom{dfrac{a}{b}} end{bmatrix},,quad a_{ij} in mathbb C tag{13}label{13} end{equation} Now we have begin{equation} rm P_{psi}boldsymbol{=}|psiranglelanglepsi|boldsymbol{=} begin{bmatrix} :z_1 : vphantom{dfrac{a}{b}}: :z_2 : vphantom{dfrac{a}{b}}: :z_3 : vphantom{dfrac{a}{b}}: end{bmatrix} begin{bmatrix} :overline{z_1} : & :overline{z_2} : & :overline{z_3} : vphantom{dfrac{a}{b}}: end{bmatrix}boldsymbol{=} begin{bmatrix} z_1overline{z_1} & z_1overline{z_2} & z_1overline{z_3} vphantom{dfrac{a}{b}}: z_2overline{z_1} & z_2overline{z_2} & z_2overline{z_3} vphantom{dfrac{a}{b}}: z_3overline{z_1} & z_3overline{z_2} & z_3overline{z_3} vphantom{dfrac{a}{b}}: end{bmatrix} tag{14}label{14} end{equation} We use the symbol $:rm P_{psi}:$ for the operator $:|psiranglelanglepsi|:$ because if $:|psirangle:$ is a normalized state vector,$:langlepsi|psirangleboldsymbol{=}1$, then $:rm P_{psi}:$ is the $''$projector$''$ on the direction $:|psirangle:$. The operator $:A,|psiranglelanglepsi|:$ is represented by the matrix begin{equation} A,|psiranglelanglepsi|boldsymbol{=} begin{bmatrix} a_{11} & a_{12} & a_{13} vphantom{dfrac{a}{b}}: a_{21} & a_{22} & a_{23} vphantom{dfrac{a}{b}}: a_{31} & a_{32} & a_{33} vphantom{dfrac{a}{b}} end{bmatrix}, begin{bmatrix} z_1overline{z_1} & z_1overline{z_2} & z_1overline{z_3} vphantom{dfrac{a}{b}}: z_2overline{z_1} & z_2overline{z_2} & z_2overline{z_3} vphantom{dfrac{a}{b}}: z_3overline{z_1} & z_3overline{z_2} & z_3overline{z_3} vphantom{dfrac{a}{b}}: end{bmatrix} tag{15}label{15} end{equation} that is ^{begin{equation} A,|psiranglelanglepsi|boldsymbol{=} begin{bmatrix} begin{array}{c|c|c} a_{11}z_1overline{z_1}boldsymbol{+}a_{12}z_2overline{z_1}boldsymbol{+}a_{13}z_3overline{z_1} & cdots & cdots vphantom{dfrac{tfrac{a}{b}}{tfrac{a}{b}}} hline cdots & a_{21}z_1overline{z_2}boldsymbol{+}a_{22}z_2overline{z_2}boldsymbol{+}a_{23}z_3overline{z_2} & cdots vphantom{dfrac{tfrac{a}{b}}{tfrac{a}{b}}} hline cdots & cdots & a_{31}z_1overline{z_3}boldsymbol{+}a_{32}z_2overline{z_3}boldsymbol{+}a_{33}z_3overline{z_3} vphantom{dfrac{tfrac{a}{b}}{tfrac{a}{b}}} end{array} end{bmatrix} tag{16}label{16} end{equation}} We show only the diagonal elements since we are interested for their sum begin{align} texttt{Tr}left(A|psiranglelanglepsi|right) boldsymbol{=} & hphantom{ ...}left(a_{11}z_1overline{z_1}boldsymbol{+}a_{12}z_2overline{z_1}boldsymbol{+}a_{13}z_3overline{z_1}right) nonumber & boldsymbol{+}left(a_{21}z_1overline{z_2}boldsymbol{+}a_{22}z_2overline{z_2}boldsymbol{+}a_{23}z_3overline{z_2}right) nonumber & boldsymbol{+}left(a_{31}z_1overline{z_3}boldsymbol{+}a_{32}z_2overline{z_3}boldsymbol{+}a_{33}z_3overline{z_3}right) nonumber boldsymbol{=} & hphantom{ ...,}overline{z_1}left(a_{11}z_1boldsymbol{+}a_{12}z_2boldsymbol{+}a_{13}z_3right) nonumber & boldsymbol{+}overline{z_2}left(a_{21}z_1boldsymbol{+}a_{22}z_2boldsymbol{+}a_{23}z_3right) nonumber & boldsymbol{+}overline{z_3}left(a_{31}z_1boldsymbol{+}a_{32}z_2boldsymbol{+}a_{33}z_3right) nonumber boldsymbol{=} & begin{bmatrix} :overline{z_1} : & :overline{z_2} : & :overline{z_3} : vphantom{dfrac{a}{b}}: end{bmatrix} begin{bmatrix} :a_{11}z_1boldsymbol{+}a_{12}z_2boldsymbol{+}a_{13}z_3 : vphantom{dfrac{a}{b}}: :a_{21}z_1boldsymbol{+}a_{22}z_2boldsymbol{+}a_{23}z_3 : vphantom{dfrac{a}{b}}: :a_{31}z_1boldsymbol{+}a_{32}z_2boldsymbol{+}a_{33}z_3 : vphantom{dfrac{a}{b}}: end{bmatrix} nonumber boldsymbol{=} & begin{bmatrix} :overline{z_1} : & :overline{z_2} : & :overline{z_3} : vphantom{dfrac{a}{b}}: end{bmatrix} begin{bmatrix} a_{11} & a_{12} & a_{13} vphantom{dfrac{a}{b}}: a_{21} & a_{22} & a_{23} vphantom{dfrac{a}{b}}: a_{31} & a_{32} & a_{33} vphantom{dfrac{a}{b}} end{bmatrix} begin{bmatrix} :z_1 : vphantom{dfrac{a}{b}}: :z_2 : vphantom{dfrac{a}{b}}: :z_3 : vphantom{dfrac{a}{b}}: end{bmatrix} nonumber boldsymbol{=} & :langlepsi|A|psirangle tag{17}label{17} end{align} qed.

$=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!$

Finally we run the proof for an infinite-dimensional separable Hilbert space $:mathbb H:$ with infinite countable complete orthonormal basis $:{phi_mu}:$ having in mind (or under the table) the finite case above.

First we must say that if $:mathcal O:$ is an operator on the Hilbert space $:mathbb H:$ then its trace is defined by begin{equation} texttt{Tr}left(mathcal O right) boldsymbol{equiv} langlephi_mu|mathcal O|phi_murangle , quad texttt{(Einstein summation convention)} tag{18}label{18} end{equation} The trace is well-defined : it actually depends on $:mathcal O:$ and not on the basis $:{phi_mu}:$ used. In other words it's invariant.

The states $:|psirangle:$ and $:langlepsi|:$ are expressed in the $:{phi_mu}:$ basis as follows
begin{align} |psirangle & boldsymbol{=} |phi_murangle:langlephi_mu|psirangle tag{19a}label{19a} langlepsi| & boldsymbol{=} langlephi_nu|:overline{langlephi_nu|psirangle} tag{19b}label{19b} end{align} corresponding to equations eqref{11b},eqref{11d} of the finite 3-dimensional case, while the complex numbers $:langlephi_mu|psirangle,overline{langlephi_nu|psirangle}:$ correspond to the coordinates $:z_mu,overline{z_nu}$.

So begin{equation} rm P_{psi}boldsymbol{=}|psiranglelanglepsi|boldsymbol{=} Bigl(|phi_muranglelanglephi_mu|psirangleBigr) Bigl(langlephi_nu|overline{langlephi_nu|psirangle}Bigr)quadboldsymbol{implies} nonumber end{equation} begin{equation} rm P_{psi}boldsymbol{=}|psiranglelanglepsi|boldsymbol{=} Bigl(langlephi_mupsirangle: overline{langlephi_nu|psirangle}vphantom{dfrac{a}{b}}Bigr) :|phi_muranglelanglephi_nu|boldsymbol{=} |phi_murangleBigl(langlephi_mupsirangle: overline{langlephi_nu|psirangle}vphantom{dfrac{a}{b}}Bigr) :langlephi_nu| tag{20}label{20} end{equation} Formally the operator $:rm P_{psi}:$ could be represented by a $''$square$''$ matrix of infinite but countable rows and columns as follows begin{equation} rm P_{psi}boldsymbol{=}|psiranglelanglepsi|boldsymbol{=} begin{bmatrix} begin{array}{c|c|c|c|c} langlephi_1|psirangle |overline{langlephi_1|psirangle} & langlephi_1|psirangle |overline{langlephi_2|psirangle} & boldsymbol{cdots} & langlephi_1|psirangle |overline{langlephi_nu|psirangle} & boldsymbol{cdots}vphantom{dfrac{a}{b}}: langlephi_2|psirangle |overline{langlephi_1|psirangle} & langlephi_2|psirangle |overline{langlephi_2|psirangle} & boldsymbol{cdots} & langlephi_2|psirangle |overline{langlephi_nu|psirangle} & boldsymbol{cdots}vphantom{dfrac{a}{b}}: boldsymbol{cdots} & boldsymbol{cdots} & boldsymbol{cdots} & boldsymbol{cdots} & boldsymbol{cdots} vphantom{dfrac{a}{b}}: langlephi_mu|psirangle |overline{langlephi_1|psirangle} & langlephi_mu|psirangle |overline{langlephi_2|psirangle} & boldsymbol{cdots} & langlephi_mu|psirangle |overline{langlephi_nu|psirangle} & boldsymbol{cdots}vphantom{dfrac{a}{b}}: boldsymbol{cdots} & boldsymbol{cdots} & boldsymbol{cdots} & boldsymbol{cdots} & boldsymbol{cdots} vphantom{dfrac{a}{b}}: end{array} end{bmatrix} tag{21}label{21} end{equation} corresponding to the matrix in the last to the right side of equation eqref{14} for the finite case.

From the definition of the trace, see equation eqref{18} begin{align} texttt{Tr}left(A|psiranglelanglepsi|right) boldsymbol{=} & Big<phi_nuvphantom{dfrac{a}{b}}Big| A|psiranglelanglepsi|Big|phi_nuBig>boldsymbol{=}Big<phi_nuvphantom{dfrac{a}{b}}Big| A,Big|psilanglepsi|phi_nurangleBig>boldsymbol{=}Big<phi_nuvphantom{dfrac{a}{b}}Big| A,Big| overline{langlephi_nu|psirangle}psiBig> nonumber boldsymbol{=} & Big<phi_nuoverline{langlephi_nu|psirangle}vphantom{dfrac{a}{b}}Big| A,Big| psiBig>boldsymbol{=}langlepsi|A|psirangle tag{22}label{22} end{align} QED.

$=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!$

begin{align} & overbrace{begin{bmatrix} boldsymbol{longleftarrow} langlephi_nu| boldsymbol{longrightarrow} vphantom{tfrac{a}{b}} end{bmatrix}boldsymbol{cdot} begin{bmatrix} boldsymbol{nwarrow} & boldsymbol{uparrow} & boldsymbol{nearrow} vphantom{dfrac{a}{b}}: boldsymbol{leftarrow} & A & boldsymbol{rightarrow} vphantom{dfrac{a}{b}}: boldsymbol{swarrow} & boldsymbol{downarrow} & boldsymbol{searrow} vphantom{dfrac{a}{b}} end{bmatrix}boldsymbol{cdot} begin{bmatrix} :boldsymbol{uparrow} : vphantom{dfrac{a}{b}}: :|psirangle : vphantom{dfrac{a}{b}}: :boldsymbol{downarrow} : vphantom{dfrac{a}{b}}: end{bmatrix}boldsymbol{cdot} underbrace{ begin{bmatrix} boldsymbol{longleftarrow} langlepsi|boldsymbol{longrightarrow} vphantom{tfrac{a}{b}} end{bmatrix}boldsymbol{cdot} begin{bmatrix} :boldsymbol{uparrow} : vphantom{dfrac{a}{b}}: :|phi_{nu}rangle : vphantom{dfrac{a}{b}}: :boldsymbol{downarrow} : vphantom{dfrac{a}{b}}: end{bmatrix}}_{texttt{scalar}boldsymbol{=}langlepsi|phi_nurangle boldsymbol{=}overline{langlephi_nu|psirangle}}}^{Big<phi_nuvphantom{dfrac{a}{b}}Big| A|psiranglelanglepsi|Big|phi_nuBig>} nonumber & tag{23}label{23} &boldsymbol{=} underbrace{underbrace{begin{bmatrix} boldsymbol{longleftarrow} langlephi_nu|,overline{langlephi_nu|psirangle}boldsymbol{longrightarrow} vphantom{tfrac{a}{b}} end{bmatrix}}_{langlepsi|} boldsymbol{cdot} begin{bmatrix} boldsymbol{nwarrow} & boldsymbol{uparrow} & boldsymbol{nearrow} vphantom{dfrac{a}{b}}: boldsymbol{leftarrow} & A & boldsymbol{rightarrow} vphantom{dfrac{a}{b}}: boldsymbol{swarrow} & boldsymbol{downarrow} & boldsymbol{searrow} vphantom{dfrac{a}{b}} end{bmatrix}boldsymbol{cdot} begin{bmatrix} :boldsymbol{uparrow} : vphantom{dfrac{a}{b}}: :|psirangle : vphantom{dfrac{a}{b}}: :boldsymbol{downarrow} : vphantom{dfrac{a}{b}}: end{bmatrix}}_{langlepsi|A|psirangle} nonumber end{align}

Note that $x=langlepsi|hat{A}|psirangle$ is just a number. (I dropped the subscript $a$ as it does not seem to play any role.) Therefore, tr${x}=x$. The trace of an operator $hat{O}$ is formally defined in terms of any complete orthogonal basis $|nrangle$ as $$ text{tr}{hat{O}} = sum_n langle n|hat{O}|nrangle . $$ Any complete basis can be used to resolve the identity $$ mathbf{1} = sum_n |nranglelangle n| . $$ Now we can insert the resolved identity into the expression for $x$: $$ x=sum_nlanglepsi|nranglelangle n|hat{A}|psirangle = sum_nlangle n|hat{A}|psiranglelanglepsi|nrangle = text{tr}{hat{A}|psiranglelanglepsi|} . $$ All we did was to exchange the two numbers under the summation to obtain the formal expression for the trace.