A real $n\times n$ symmetric matrix $X=(x_{ij})_{n \times n}$ is
called a bisymmetric matrix if $x_{ij}=x_{n+1-j, n+1-i}$. Based on
the projection theorem, the canonical correlation decomposition and
the generalized singular value decomposition, a method useful for
finding the least-squares solutions of the matrix equation
$A^{T}XA=B$ over bisymmetric matrices is proposed. The expression of
the least-squares solutions is given. Moreover, in the corresponding
solution set, the optimal approximate solution to a given matrix is
also derived. A numerical algorithm for finding the optimal
approximate solution is also described.