A TikZ based solution:
\documentclass{article}
\usepackage{amsmath,amssymb}
\usepackage{tikz}
\usetikzlibrary{matrix,calc}
\begin{document}
\begin{align*}
P&=
\begin{tikzpicture}[baseline=-\the\dimexpr\fontdimen22\textfont2\relax ]
\matrix(m)[matrix of math nodes,left delimiter=(,right delimiter=),inner sep=4pt,ampersand replacement=\&]
{
x_1 \& y_1 \& s_1 \& z_1 \\
x_2 \& y_2 \& s_2 \& z_2 \\
x_3 \& y_3 \& s_3 \& z_3 \\
x_4 \& y_4 \& s_4 \& z_4 \\
};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\foreach \s/\stext in {m-4-1/1,m-4-2/2,m-4-3/3,m-4-4/4}{
% bottom index
\node[red,shift=(\s.south),yshift=-0.4cm](0,0) {$a_{\stext}$} ;
}
\foreach \n/\ntext in {m-1-1/1,m-1-2/2,m-1-3/3,m-1-4/4}{
% top index
\node[blue,shift=(\n.north),yshift=0.4cm](0,0) {$b_{\ntext}$} ;
}
\end{tikzpicture}
\end{align*}
\end{document}

IMPROVEMENT
In his comment percusse suggested a very nice way of make foreachs simpler: here is the previous example enhanced.
\documentclass{article}
\usepackage{amsmath,amssymb}
\usepackage{tikz}
\usetikzlibrary{matrix,calc}
\begin{document}
\begin{align*}
P&=
\begin{tikzpicture}[baseline=-\the\dimexpr\fontdimen22\textfont2\relax ]
\matrix(m)[matrix of math nodes,left delimiter=(,right delimiter=),inner sep=4pt,ampersand replacement=\&]
{
x_1 \& y_1 \& s_1 \& z_1 \\
x_2 \& y_2 \& s_2 \& z_2 \\
x_3 \& y_3 \& s_3 \& z_3 \\
x_4 \& y_4 \& s_4 \& z_4 \\
};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\foreach \s in {1,2,...,4}{
% bottom index
\node[red,shift=(m-4-\s.south),yshift=-0.4cm](0,0) {$a_{\s}$};
}
\foreach \n in {1,2,...,4}{
% top index
\node[blue,shift=(m-1-\n.north),yshift=0.4cm](0,0) {$b_{\n}$} ;
}
\end{tikzpicture}
\end{align*}
\end{document}
blkarray; there are examples on this site. – egreg Jun 18 '12 at 13:38