The real answer to this question is that if you play with it long enough, you'll hit an algorithm requiring 7 multiplications – not necessarily the same as Strassen's, but an equivalent one, in a certain sense: it is known that all such algorithms are equivalent, as shown by de Groote in his 1978 paper, On varieties of optimal algorithms for the computation of bilinear mappings. II. Optimal algorithms for 2 × 2-matrix multiplication.
There are many attempts in the literature to explain how one could come up with such an algorithm, for example:
- Gideon Yuval, A simple proof of Strassen’s result, 1978. We explain this approach below.
- Ann Q. Gates, Vladik Kreinovich, Strassen's Algorithm Made (Somewhat) More Natural: A Pedagogical Remark, 2001. The idea is to use symmetries to guess the linear combinations corresponding to one of the matrices being multiplied, and then to pair them intelligently with linear combinations of the other matrix.
- Jacob Minz, Derivation of Strassen's Algorithm for the multiplication of 2×2 matrices, 2015. The idea is to apply linear transformations to obtain a simpler looking problem, and then to solve it by hand.
- Christian Ikenmeyer, Vladimir Lysikov, Strassen's 2x2 matrix multiplication algorithm: A conceptual perspective, 2017. The idea is to consider a basis for the $2\times 2$ traceless matrices, and use its multiplication table to construct Strassen's algorithm.
- Joshua A. Grochow, Christopher Moore, Designing Strassen’s Algorithm, 2017. The idea is to start with vectors which form the vertices of an equilateral triangle in the plane, and use elementary properties of these vectors to come up with an algorithm.
Gideon Yuval shows how you could come up with Strassen's algorithm. The starting point is to convert matrix multiplication to the problem of computing a matrix-vector product: computing
$$ \begin{pmatrix} a & b \\ c & d \end{pmatrix} \times \begin{pmatrix} e & f \\ g & h \end{pmatrix} $$
is the same as computing
$$
\begin{pmatrix}
a & 0 & b & 0 \\
0 & a & 0 & b \\
c & 0 & c & 0 \\
0 & d & 0 & d
\end{pmatrix}
\times
\begin{pmatrix} e \\ f \\ g \\ h \end{pmatrix}
$$
Suppose that we could write the matrix on the left as a sum $\ell_1 M_1 + \cdots + \ell_7 M_7$, where $\ell_i$ is a linear combination of $a,b,c,d$ and $M_i$ is a rank one matrix, say $M_i = x_i y_i^T$. The product we are after is thus
$$
\sum_{i=1}^7 \ell_i M_i \begin{pmatrix} c\\d\\e\\f \end{pmatrix} =
\sum_{i=1}^7 \ell_i x_i y_i^T \begin{pmatrix} c\\d\\e\\f \end{pmatrix} =
\sum_{i=1}^7 \ell_i r_i x_i,
$$
where $r_i$ is a linear combination of $e,f,g,h$. This shows that each entry of the product matrix is some linear combination of the products $\ell_i,r_i$.
Let us now show how one could find the decomposition. We start by cancelling the top left and bottom right entries, in a way which avoids hitting zero entries:
$$
\begin{pmatrix}
a & 0 & b & 0 \\
0 & a & 0 & b \\
c & 0 & d & 0 \\
0 & c & 0 & d
\end{pmatrix} =
\begin{pmatrix}
a & 0 & a & 0 \\
0 & 0 & 0 & 0 \\
a & 0 & a & 0 \\
0 & 0 & 0 & 0
\end{pmatrix} +
\begin{pmatrix}
0 & 0 & 0 & 0 \\
0 & d & 0 & d \\
0 & 0 & 0 & 0 \\
0 & d & 0 & d
\end{pmatrix} +
\begin{pmatrix}
0 & 0 & b-a & 0 \\
0 & a-d & 0 & b-d \\
c-a & 0 & d-a & 0 \\
0 & c-d & 0 & 0
\end{pmatrix}
$$
This results in a mess, which we try to fix by "flipping" the inner square:
$$
\begin{pmatrix}
0 & 0 & b-a & 0 \\
0 & a-d & 0 & b-d \\
c-a & 0 & d-a & 0 \\
0 & c-d & 0 & 0
\end{pmatrix} =
\begin{pmatrix}
0 & 0 & 0 & 0 \\
0 & a-d & a-d & 0 \\
0 & d-a & d-a & 0 \\
0 & 0 & 0 & 0
\end{pmatrix} +
\begin{pmatrix}
0 & 0 & b-a & 0 \\
0 & 0 & d-a & b-d \\
c-a & a-d & 0 & 0 \\
0 & c-d & 0 & 0
\end{pmatrix}
$$
Since $d-a = (b-a)-(b-d)$ and $a-d = (c-d)-(c-a)$, it is easy to represent the last matrix as a sum of four rank one matrices:
$$
\begin{pmatrix}
0 & 0 & b-a & 0 \\
0 & 0 & d-a & b-d \\
c-a & a-d & 0 & 0 \\
0 & c-d & 0 & 0
\end{pmatrix} =
\begin{pmatrix}
0 & 0 & b-a & 0 \\
0 & 0 & b-a & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0
\end{pmatrix} +
\begin{pmatrix}
0 & 0 & 0 & 0 \\
0 & 0 & d-b & b-d \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0
\end{pmatrix} +
\begin{pmatrix}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & c-d & 0 & 0 \\
0 & c-d & 0 & 0
\end{pmatrix} +
\begin{pmatrix}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
c-a & a-c & 0 & 0 \\
0 & 0 & 0 & 0
\end{pmatrix}
$$
In total, we obtain the following representation:
$$
\begin{pmatrix}
a & 0 & b & 0 \\
0 & a & 0 & b \\
c & 0 & d & 0 \\
0 & c & 0 & d
\end{pmatrix} =
\begin{pmatrix}
a & 0 & a & 0 \\
0 & 0 & 0 & 0 \\
a & 0 & a & 0 \\
0 & 0 & 0 & 0
\end{pmatrix} +
\begin{pmatrix}
0 & 0 & 0 & 0 \\
0 & d & 0 & d \\
0 & 0 & 0 & 0 \\
0 & d & 0 & d
\end{pmatrix} +
\begin{pmatrix}
0 & 0 & 0 & 0 \\
0 & a-d & a-d & 0 \\
0 & d-a & d-a & 0 \\
0 & 0 & 0 & 0
\end{pmatrix} +
\begin{pmatrix}
0 & 0 & b-a & 0 \\
0 & 0 & b-a & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0
\end{pmatrix} +
\begin{pmatrix}
0 & 0 & 0 & 0 \\
0 & 0 & d-b & b-d \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0
\end{pmatrix} +
\begin{pmatrix}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & c-d & 0 & 0 \\
0 & c-d & 0 & 0
\end{pmatrix} +
\begin{pmatrix}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
c-a & a-c & 0 & 0 \\
0 & 0 & 0 & 0
\end{pmatrix}
$$
