  We are presenting an Internal Linear Combination (ILC) CMB map, in which the
foreground is reduced through harmonic variance minimization. We have derived
our method by converting a general form of pixel-space approach into spherical
harmonic space, maintaining full correspondence. By working in spherical
harmonic space, spatial variability of linear weights is incorporated in a
self-contained manner and our linear weights are continuous functions of
position over the entire sky. The full correspondence to pixel-space approach
enables straightforward physical interpretation on our approach. In variance
minimization of a linear combination map, the existence of a cross term between
residual foregrounds and CMB makes the linear combination of minimum variance
differ from that of minimum foreground. We have developed an iterative
foreground reduction method, where perturbative correction is made for the
cross term. Our CMB map derived from the WMAP data is in better agreement with
the WMAP best-fit $\Lambda$CDM model than the WMAP team's Internal Linear
Combination map. We find that our method's capacity to clean foreground is
limited by the availability of enough spherical harmonic coefficients of good
Signal-to-Noise Ratio (SNR).
