  We analyze the $B \to \phi K^{*}$ polarization puzzle in the Minimal
Supersymmetric Standard Model (MSSM) including the neutral Higgs boson (NHB)
contributions. To calculate the non-factorizable contributions to hadronic
matrix elements of operators, we have used the QCD factorization framework to
the $\alpha_s $ order. It is shown that the recent experimental results of the
polarization fractions in $B\to \phi K^{*}$ decays, which are difficult to be
explained in SM, could be explained in MSSM if there are flavor non-diagonal
squark mass matrix elements of 2nd and 3rd generations, which also satisfy all
relevant constraints from known experiments ($B\to X_s\gamma, B_s\to
\mu^+\mu^-, B\to X_s \mu^+\mu^-, B\to X_s g, \Delta M_s$, etc.). We have shown
in details that the experimental results can be accommodated with the flavor
non-diagonal mass insertion of chirality RL, RL+LR, RR, or LL+ RR when the NHB
contributions as well as $\mathcal{O}(\alpha_s)$ corrections of hadronic matrix
elements of operators are included. However the branching ratios for the decay
are smaller than the experimental measurements.
