Summary:
A dimension restriction is necessary. Lemma 2 below observes that if arbitrary dimension is allowed the problem (even restricted to permutation matrices) is at least as hard as Graph Isomorphism (GI-hard).
OP's problem (as restricted by OP to instances having constant dimension and matrices that either (i) all have at most one 1 per row, or (ii) all have at most one 1 per column) has a poly-time algorithm.
The upper bound we show is actually slightly stronger.
Given an instance of OP's problem,
let $\mathcal A$ denote the set of matrices labeling the edges.
Let $\mathcal A^*\vec c$ be the set of vectors obtainable
as products of the form $M_1 M_2 \cdots M_k \vec c$ for any $k\ge 0$, with each $M_i\in \mathcal A$.
Lemma 1. The problem has a poly($|G|+|\mathcal A^*\vec c|)$-time algorithm (with no restriction on the matrices or dimension).
For example, the following special cases are in P:
If every edge has the identity matrix $I$, then $\mathcal A^*\vec c=\{\vec c\}$, so $|\mathcal A^*\vec c| = 1$. (This is the result shown in D.W.'s answer.)
If the matrices are (0,1) matrices, either all having at most one 1 in each column, or all having at most one 1 in each row,
then $|\mathcal A^*\vec c|\le |\mathcal A^*| \le n^n$.
So when the dimension $n$ is constant (or even $O(\log(|G|)/\log\log |G|)$) the algorithm
runs in time poly$(|G|)$.
Proof sketch for Lemma 1.
Given the digraph $G=(V,E)$,
the algorithm computes $\mathcal A^*\vec c$.
It then constructs a bigger edge-weighted digraph $G'$
with a vertex $q(v, y)$
for every pair $(v, y)$ in $V \times \mathcal A^*\vec c$.
For every pair of vertices $q(v, y)$ and $q(v', y')$,
there is an edge from the former to the latter
iff
- $(v, v')$ is in $E$, and
- $y = Ay'$, where $A$ is the matrix labeling edge $(v, v')$.
If there is such an edge,
it is given (scalar) weight $\vec b\, \vec y'$.
The algorithm then uses any standard shortest-path algorithm
to find the shortest path from any starting vertex in $G'$
to the ending vertex in $G'$,
where the starting vertices in $G'$
are those of the form $q(u, y)$
where $u$ is the starting vertex for OP's instance
and $y\in\mathcal A^*\vec c$,
and the ending vertex is $q(v, \vec c)$,
where $v$ is the ending vertex for OP's instance.
The idea for correctness is that any path
$q(v_1, y_1), q(v_2, y_2), \ldots, q(v_k, y_k=\vec c)$
in $G'$
from a starting vertex to the ending vertex
represents the path
$(v_1, v_2, \ldots, v_k)$ in $G$
from the starting vertex $v_1$ to the ending vertex $v_k$ in $G$.
The conditions for edges in $G'$
guarantee that after this path traverses an edge
$(v_i, v_{i+1})$ labeled with some matrix $A$ and vector $\vec b$,
it is the case that the product
$M_{i+1} M_{i+2} \cdots M_k \vec c$
along the remainder of the path will equal $y_{i+1}$.
Hence, (in $G)$ the vector $\vec b$ of the edge $(v_i, v_{i+1})$
will contribute $\vec b y_{i+1}$ to the final cost,
just as the corresponding edge
(with weight $\vec b y_{i+1}$) on the path in $G'$
contributes to that path's final cost.
Thus, the cost of the path in $G'$
will equal the cost of the corresponding path in $G$.
$~~~~\Box$
Lemma 2. The problem (restricted to permutation matrices,
but not constant dimension) is GI-hard.
Proof sketch for Lemma 2. The proof is by a reduction from Graph Isomorphism (GI).
Given a GI instance $(G_1, G_2)$ of two undirected graphs, each (WLOG) with vertex set $[n]$, and with the same number of edges, the reduction outputs the instance of OP's problem defined as follows. (Note that here $n$ is the number of vertices in $G_1$
and $G_2$, not the dimension of an instance of OP's problem.
The instance of OP's problem produced by the reduction
will have dimension $n\choose 2$.)
Let $U$ denote the set of possible edges
(i.e. unordered pairs of distinct vertices)
from vertex set $[n]$.
So $|U| = {n\choose 2}$ and, for example,
the edge sets of $G_1$ and $G_2$
are subsets of $U$.
Let $\{0,1\}^U$
represent the indicator vectors for such subsets of $U$
in the standard way:
a given $x\in\{0,1\}^U$
is identified with the edge set $E$ such that $x_e = 1$ iff $e\in E$.
For $(u,v)\in U$,
let $\pi(u, v)$ denote the permutation
induced on $U$ by exchanging the (names of the) vertices $u$ and $v$.
That is, for every vertex $w\not\in \{u, v\}$,
$\pi(u,v)$ swaps $(u, w) \in U$ with $(v, w)\in U$,
leaving all other possible edges unchanged.
Let $\Pi = \{\pi(u, v) : (u,v)\in U\}$ denote the set of such permutations. Note $|\Pi| = {n\choose 2}$.
Now $G_1=([n], E_1)$ is isomorphic to $G_2=([n], E_2)$
if and only if there is a sequence of
permutations in $\Pi$
such that applying the sequence of permutations
to (the indicator vector $x(E_1) \in \{0,1\}^U$ of) $E_1$
yields (the indicator vector $x(E_2) \in \{0,1\}^U$ of) $E_2$.
Furthermore, there is a constant $c>0$
such that there is such a sequence
if and only if there is such a sequence of length
at most $cn \log n$.
The reduction outputs a multi-digraph $D$
with vertex set $[N]$ where $N=\lceil c n\log n\rceil$
and, for each $i\in [N-1]$,
there are $|\Pi|$ multi-edges from $i$ to $i+1$,
with each such edge being labeled with the matrix for a (distinct) permutation in $\Pi$.
Thus, the paths from $1$ to $N$ in $D$ correspond
to the sequences of permutations in $\Pi$.
The start vertex $1$ is labeled
with (the starting) vector $x(E_1)$.
The target vertex $N$ is labeled
with the vector $\vec c = -x(E_2)$.
Then, the cost of taking a path from $1$ to $N$
is (at most) $-|E_1| = -|E_2|$ iff
the sequence of permutations induced by the path,
applied to $x(E_1)$, yields $x(E_2)$.
Thus, this cost is achievable
if and only if $G_1$ and $G_2$ are isomorphic.
(Note that if desired the multi-digraph $D$ can be converted
into an equivalent digraph using the standard
idea of replacing each edge $(i, j)$
by a path $(i, v, j)$,
where $v$ is a new artificial vertex, unique to the edge.)
$~~~\Box$
