This is a typical deconvolution problem that you can solve either by transforming to the frequency domain where convolution is a simple multiplication:
\begin{equation}
F\left \{ {\bf{I}} \right \}= F\left \{ {\bf{z}} \right \} F\left \{ {\bf{f}} \right \},
\end{equation}
where ${\mathcal{F}}$ denotes the Fourier transform (or DFT) so $f$ will be:
\begin{equation}
{\bf{f}}={\mathcal{F}}^{-1}\left \{ \frac{{\mathcal{F}}\left \{ {\bf{I}} \right \}}{{\mathcal{F}}\left \{ {\bf{z}} \right \}} \right \}
\end{equation},
where ${\mathcal{F}}^{-1}$ is the inverse Fourier transform.
Otherwise, you can also represent in a matrix form, using a Circulant (kind of Toeplitz) for the convolution operation. The frequency domain calculation though does not involve inverting (or pseudo inverting) a matrix.
If there's additive noise in your problem, please see Wiener deconvolution (both time and frequency domain implementation)