  In this paper, we develop a method for solving the problem of minimizing the
$H^2$ error norm between the transfer functions of original and reduced systems
on the set of stable matrices and two Euclidean spaces. That is, we develop a
method for identifying the optimal reduced system from all stable linear
systems. However, it is difficult to develop an algorithm for solving this
problem, because the set of stable matrices is highly non-convex. To overcome
this issue, we show that the problem can be transformed into a tractable
Riemannian optimization on the product manifold of the set of skew-symmetric
matrices, the manifold of the symmetric positive-definite matrices, and two
Euclidean spaces. The stability of the reduced systems constructed using the
optimal solutions to our problem is preserved. To solve the reduced problem,
the Riemannian gradient and Hessian are derived and a Riemannian trust-region
method is developed. The initial point in the proposed approach is selected
using the output from the balanced truncation (BT) method. Numerical
experiments demonstrate that our method considerably improves the results given
by BT in the sense of the $H^2$ norm, and also provides reduced systems that
are globally near-optimal solutions to the problem of minimizing the
$H^{\infty}$ error norm. Moreover, we show that our method provides a better
reduced model than BT from the viewpoint of the frequency response.
