I am interesting in some sort of algorithm for calculating the Boolean rank of small $M \times N$ Boolean matrices.
Just to be clear, by Boolean matrices I mean matrices with entries $0$ or $1$ where $1+1=1$, and I am defining the Boolean rank as:
Matrix $X$ has $rank(X) = 1$ iff $X = ab^T$ for column vectors $a$ and $b$
Matrix $X$ has $rank(X) ≤ k$ if it can be represented as a sum of $k$ rank-1 matrices. Smallest such $k$ is the rank of $X$
For example,
$A = \pmatrix{1&1&0\\1&1&1\\0&1&1} = \pmatrix{1\\1\\0}\pmatrix{1&1&0} + \pmatrix{0\\1\\1}\pmatrix{0&1&1} = \pmatrix{1&1&0\\1&1&0\\0&0&0}+\pmatrix{0&0&0\\0&1&1\\0&1&1}$
so $rank(A)=2$
I am also interested in some way to factor an $M \times N$ Boolean matrix with rank $k$ into $M \times k$ and $k \times N$ Boolean matrices.
For example,
$A = \pmatrix{1&1&0\\1&1&1\\0&1&1} =\pmatrix{1&0\\1&1\\0&1}\pmatrix{1&1&0\\0&1&1}$
For a matrix $X$ of rank $r$ I will call the right hand side of $$X = a_1b_1^T + \dots + a_rb_r^T$$ a $decomposition$ of $X$. I will call $a_ib_i^T$ a term in the decomposition.
Any rank 1-matrix will be a matrix where a submatrix is made entirely out of ones and the rest of the entries are zeros. A submatrix $B'$ of the matrix $B$ is a matrix where we pick out some columns and some rows. E.g. for the matrix $$B = \begin{pmatrix} 1 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \\ 1 & 1 & 1 & 0 \end{pmatrix}$$ the submatrix $B'$ consisting of rows 1 and 3 and columns 2 and 3 is: $$B ' = \begin{pmatrix} 0 & 0 \\ 1 & 1 \end{pmatrix}.$$
You can reformulate your problem to other problems:
Reformulate to a set problem
This can be reformulated to a set problem. Given an $M \times N$ binary matrix $B$, your universe $U$ is the set of indices $(1,1), (1, 2), \dots, (M,N)$, and you have a set $\mathfrak B \subseteq U$ consisting of the indices where $B$ is nonzero: $$\mathfrak B = \{ (i,j) : B_{ij} = 1 \}$$ You have a collection of sets $\mathcal R_1$, where each $R \in \mathcal R_1$ corresponds to a rank 1 matrix, the set $R$ consisting of the indices where the matrix is nonzero.
So, then your problem is to find a smallest $\mathcal C \subset R$ such that the union of all sets in $\mathcal C$ is equal to $\mathfrak B$, i.e. the $\mathcal C$ with smallest size $| C |$ such that $$\bigcup_{C \in \mathcal C} C = \mathfrak B.$$ This, I think, is a special case of the set basis problem.
Reformulate to a graph problem
You can also, given an $M \times N$ binary matrix $B$, construct a bipartite graph $G_B$ consisting of two vertex sets $V_1, V_2$ of size $M, N$ respectively, with an edge between $v_i \in V_1$ and $w_j \in V_2$ if $B_{ij} = 1$, i.e. the elements in $V_1$ correspond to the row numbers and the elements in $V_2$ correspond to the column numbers.
The rank is then the least number of complete bipartite subgraphs needed to cover $G_B$, something that is known as the bipartite dimension of $G_B$.
You can see this by considering what a complete bipartite subgraph corresponds to: Let's say the the complete bipartite subgraph consists of the vertex sets $A, B$. Since the subgraph is complete bipartite, there is an edge between every $a \in A$ and every $b \in B$. This corresponds to a submatrix filled with ones.
Comments
Both the decision problem for the set formulation and the decision problem for the graph formulation are NP-complete.
A closely related problem, how to get as close to a given binary matrix $B$ as possible with a restricted number of terms in your decomposition, is studied in the article The Discrete Basis Problem, by Miettinen, Mielikäinen, Gionis, Das and Mannila.
An upper bound for the rank
If $B$ is an $M \times N$ boolean matrix, we see that we will never need more than $mn$ terms to decompose it. We can see this by using the standard basis. Let $e_i^n$ be the $n$-dimensional vector with a one in position $i$, zeros elsewhere. Then $e_i^n (e_j^m)^T$ is the $n \times m$ boolean matrix at position $(i,j)$. Thus we can always decompose $B$ as $$B = \sum_{i = 1}^N \sum_{j = 1}^M \alpha_{ij} e_i^N (e_j^M)^T$$ where $\alpha_{ij}$ is zero or one. We can simplify this by writing: $$B = \sum_{j = 1}^M \left( \sum_{i=1}^N \alpha_{ij} e_i^N \right)(e_j^m)^T$$ which is $B$ expressed as a sum of $M$ rank one matrices. Note that we can just as well have the $i$-sum as the outer sum, so we get that $$\operatorname{rank}(B) \leq \min \{N, M\}.$$