Algorithm to solve $L_1$ optimization of $\sum_i ||\mathbf{A_i x} - \mathbf{b_i}||_1$
I have found the answer for this problem. One way is, as suggested in the comment, to stack the matrix $\mathbf{A_i}$ and $\mathbf{b_i}$ together and apply first-order convex optimization algorithm, such as Primal-Dual Hybrid Gradient (PDHG), inexact ADMM, or their variants.
For PDHG, the problem to solve is
\begin{align}
\mathrm{min}_{\mathbf{x},\mathbf{y}} &amp;~ F(\mathbf{y}) + G(\mathbf{x}) \\
\mathrm{s.t.} &amp; ~ \mathbf{Ax} = \mathbf{y}
\end{align}
where
\begin{align}
F(\mathbf{y}) &amp;= ||\mathbf{y} - \mathbf{b}||_1 \\
G(\mathbf{x}) &amp;= 0 \\
\mathbf{A} &amp;= [\mathbf{A_1}^T, \mathbf{A_2}^T, ..., \mathbf{A_n}^T]^T \\
\mathbf{b} &amp;= [\mathbf{b_1}^T, \mathbf{b_2}^T, ..., \mathbf{b_n}^T]^T. \\
\end{align}
Thus, the convex conjugate of $F(\mathbf{y})$ is
\begin{align}
F^*(\mathbf{y}^*) &amp;= \mathrm{sup}_{\mathbf{y}} \left[\left&lt;\mathbf{y},\mathbf{y}^*\right&gt; - F(\mathbf{y})\right] \\
F^*(\mathbf{y}^*) &amp;= \left&lt;\mathbf{b},\mathbf{y^*}\right&gt; + \mathcal{I}_{\mathbf{y^*}&lt;\mathbf{b}}
\end{align}
where $\mathcal{I}_\mathcal{C}$ is an indicator function with value $0$ if the condition $\mathcal{C}$ fulfilled, otherwise $\infty$.
As the algorithm uses proximal operator, we need to get the proximal operator for $F$ and $G$,
\begin{align}
\mathrm{prox}_{F^*\sigma}(\mathbf{y}^*) &amp;= \arg\min_{\mathbf{y}'}\left\{\frac{||\mathbf{y}'-\mathbf{y}^*||^2}{2\sigma} + F^*(\mathbf{y}^*)\right\} \\
\mathrm{prox}_{F^*\sigma}(\mathbf{y}^*) &amp;= \max(\mathbf{y}^* + \sigma \mathbf{b}, \mathbf{b}) \\
\mathrm{prox}_{G\tau}(\mathbf{x}) &amp;= \mathbf{x} \\
\end{align}
Thus, the algorithm in this case reads
\begin{align}
\mathbf{y}^{n+1} &amp;= \max(\mathbf{y}^n + \sigma \mathbf{A} \mathbf{w}^n + \sigma \mathbf{b}, \mathbf{b}) \\
\mathbf{x}^{n+1} &amp;= \mathbf{x}^n - \tau \mathbf{A}^T \mathbf{y}^{n+1} \\
\mathbf{w}^{n+1} &amp;= \mathbf{x}^{n+1} + \theta (\mathbf{x}^{n+1} - \mathbf{x}^n),
\end{align}
with $\theta = 1$, and $\sigma \tau ||A||^2 &lt; 1$.
All the operations do not need the direct access to the element of the matrix $\mathbf{A}$ and only require the operation of $\mathbf{A} \cdot$ and $\mathbf{A}^T \cdot$.