This talk discusses the so-called bi-quadratic optimization over unit spheres. We show that this problem is NP-hard and there is no polynomial time algorithm returning apositive relative approximation bound. After that, we present various approximation methods based on semidefinite programming (SDP) relaxations. Our theoretical results are: For general bi-quadratic forms, we develop a $1/\max{\{m,n\}}$-approximation algorithm; for bi-quadratic forms tha are square-free, we get a relative approximation bound $1/mn$. When $\min{\{m,n\}}$ is a constant, we present two polynomial time approximation schemes (PTASs) which are based on sum of squares (SOS) relaxation hierarchy and grid sampling of the standard simplex. For practical computational purposes, we propose the first order SOS relaxation, a convex quadratic SDP relaxation and a simple minimum eigenvalue method, and give their quality analyses. Some illustrative numerical examples are also given.