Deriving the least-squares solution for this problem Suppose I have 5 images with 4000 pixels each, let this be $I$ of dimensions 4000 x 5. Suppose there exists a relationship $I = a*L$ and I want to find what $L$ is. In this case $a$ is known and has size 4000 x 3. $L$ is unknown and has size $3x5$. I want to fit the data and find the best $L$.
So first, I formulate the error function:
$$\sum_{j=1}^{4000}(I_j^T-L^Ta_j^T)^2$$
Here $I_j^T$ is a 5x1 dimension vector, think of it as a vector of 5 pixel values with each value representing one of the 5 images. $L^T$ is now size 5x3. $a_j^T$ is now size 3x1. 
Next, I take the derivative with respect to $L$ and get:
$$\frac{dE}{dL} = 2\sum_{j=1}^{4000}\left((I_j^T-L^Ta_j^T)(-a^T)\right)$$
At this point, I am lost at what to do because $(I_j^T-L^Ta_j^T)$ is size 5x1 and I can't just remove the summation by transposing it. 
Does anyone know what to do here? 
 A: Start with a clean problem statement using the familiar $\mathbf{A}x=b$ notation.
Problem Statement
Inputs
The data is given as a sequence of five images, $i_{1}$, $i_{2}$, $i_{3}$, $i_{4}$, $i_{5}$, where $i_{k}\in\mathbb{R}^{\lambda\times 1}$, with $\lambda = 4000$, for $k=1..5$.
Assume the images are linearly independent. The images must be distinct, and no image can be expressed as a combination of the other images.
Use these vectors to compose the columns of the data matrix $b$:
$$
b = \left[ \begin{array}{ccccc} i_{1} & i_{2} & i_{3} & i_{4} & i_{5} \end{array} \right] \in \mathbb{R}^{\lambda\times 5}
$$
We are given a system matrix $\mathbf{A}\in \mathbb{R}^{\lambda \times 3}$.
Output
Find the solution matrix $x_{LS}$ defined as
$$
x_{LS} = \left\{x\in \mathbb{R}^{3 \times 5}\colon \lVert \mathbf{A}x-b\rVert_{2}^{2} \text{is minimized} \right\}
$$
Verify conformability
$$
 \mathbf{A}x = b \quad \Rightarrow \quad 
\left(\lambda \times 3 \right)
\left(3 \times 5 \right)
=
\left(\lambda \times 5 \right)
$$
To connect to your notation use
$$
a \to \mathbf{A}, \quad b \to I, \quad x \to L.
$$
Solution strategies
There are many ways to find the least squares solution to $\mathbf{A}x = b$.
Ill-conditioned systems
The most powerful, yet slowest, method is the singular value decomposition.
Why does SVD provide the least squares and least norm solution to =?

How does the SVD solve the least squares problem?
A less expensive method is the $\mathbf{Q}\mathbf{R}$ composition.
Well-conditioned systems
As noted by @greg, the normal equations method may be used. But you system must be well-conditioned to extract a meaningful answer.
Deriving the Normal Equation used to solve Least Squares Problems
Deriving the least square system equations from calculus and the normal equations
Because there is no matrix $x$ such that $\mathbf{A}x-b=\mathbf{0}$, we can remedy the malady by solving an equivalent system
$$
\mathbf{A}^{T}\mathbf{A}x=\mathbf{A}^{T}b.
$$
These are the normal equations. If the $3\times3$ matrix $\left( \mathbf{A}^{T}\mathbf{A}\right)^{-1}$ exists, then the solution matrix is given by
$$
x_{LS} = \left( \mathbf{A}^{T}\mathbf{A}\right)^{-1} \mathbf{A}^{T}b.
$$
