# How do you compute the derivative of a matrix algebra expression?

I was working along with http://jimherold.com/2012/04/20/least-squares-bezier-fit/ to see if I could understand each step of the way for fitting a Bezier curve to a set of coordinates, and I understand every step except one: computing the derivative of a particular matrix multiplication.

The article gives the least squares error function between the real coordinates and the coordinates yielded by the to-be-derived Bezier function as:

$E(C_y) = \sum^n_{i=1} \left ( y_i - B(t_i) \right )^2$

Which makes a lot of sense, and then restates this as a matrix expression:

$E(C_y) = \left ( y - \mathbb{T}MC_y \right )^T \left ( y - \mathbb{T}MC_y \right )$

In this, $y$ is the y-coordinate vector, $M$ is the Bernstein polynomial coefficients for the Bezier curve in matrix form, $C_y$ is the "y coordinates we're hoping to find, to use in constructing a well-fitting curve", and $\mathbb{T}$ is a matrix of "powers of $t$ specific to fitting a curve to the polygon described by the original coordinates. Going through the article, all of that makes sense.

However, the article then describes that to find the smallest error, we want to find the roots for the derivative of the error function, and states this as solving:

$\frac{\partial E}{\partial C} = 0 = -2𝕋^T \left ( y - 𝕋MC_y \right )$

While I can verify that this works, I don't know why it works. I have no idea which rules to use to go from the original error function to this derivative. How can I get from the former to the latter?

Note that $E(C_y)$ is scalar function that depends on $c_i, i \in \{1,...n\}$ where $c_i$ are the components of $C_y$. The derivative that is being computed here is just a stacking of the partial derivative with respect to each of the $c_i's$. So, you can always expand out the expression instead of writing it in matrix-vector notation. I'll assume you don't want to do so. I'll give an answer that does indeed compute the partial for each component, but after a bit of massaging of the expression for $E(C_y)$.

Then, for this particular case (and often enough!), you will notice that $E(C_y)$ is an inner product -- I'll write this as $\langle y - TMC_y, y - TMC_y \rangle$. Over the reals, the dot product is symmetric and bilinear, and you can write this as $\langle y, y\rangle + \langle TMC_y, TMC_y \rangle - 2\langle y, TMC_y \rangle$. Each term (taken with the sign) is now a scalar function, and $E(C_y)$ is just the sum of each term. So you can split the derivative computation and add the contribution of each of the terms together as the derivative is linear ($\frac{d(f + g)}{dy} = \frac{df}{dy} + \frac{dg}{dy}$) as is the partial derivative.

At this point its best to compute for a particular component, notice the structured form that results and remember it as a 'rule'. Let's do this component wise computation for $c_1$.

The first term evaluates to a zero partial derivative for $c_1$ as it is independent of $c_1$.

Consider $x, y \in \mathbb{R}^n$ with components $x_i, y_i$ repectively, and for $A \in \mathbb{R}^{n\times n}$. Note that $\langle x, Ay \rangle = \langle A^Tx, y \rangle = y^TA^Tx = x^TAy$ as this is just a scalar, and that $\frac{\partial \langle y, x \rangle}{\partial x_1} = y_1$ (try writing this out in expanded form). Thus, the third term of the summation yields $-2 \frac{\partial \langle (TM)^Ty, C_y \rangle}{\partial c_1} = -2((TM)^Ty)_1$.

A final rule that you must work out is that $\frac{\partial \langle Ax, x \rangle}{\partial x_1} = ((A + A^T)x)_1$. I couldn't pen down a particularly nice vector-notation way of going about this, so, I leave this to you as an exercise :). When $A = A^T$, the above "rule" then reads as $2(Ax)_1$. $A^TA$ is always equal to $(A^TA)^T$. As the second term of the summation of $E(C_y)$ is $\langle (TM)^TTMC_y,C_y \rangle$, the derivative evaluates to $((TM)^TTMC_y)_1$.

Stacking up the partial derivatives, we see $\frac{\partial E}{\partial C} = -2(TM)^Ty + 2(TM)^TTMC_y$, setting which to zero, you get the so-called normal equation: $min_{C_y} ||TMC_y - y||_2 \Rightarrow C_y~\text{satisfies}~(TM)^TTMC_y = (TM)^Ty \Rightarrow T^TTMC_y = T^Ty$ as $M$ has full-rank.

Wrt. the article itself, unless I am very much mistaken, there is no reason to deem that computing a chord-length parameterized (the thing the author does when deciding the values of $t_i$) cubic Bezier curve whose control points minimize the summed square distance to the input data samples will yield a "best-fit" Bezier curve to the data, as the formulation says nothing about the behavior of the Bezier curve between the samples. This is perhaps why the author of the article didn't find "best-fit" solutions that were deterministic -- the least squares approximation is perhaps a good first cut, but not really the (subjectively defined) "best-fit".

• Absolutely, it's a "best fit" only with respect to least square fitting, after picking reasonable t values based on distance along the polygon. Those t values can be far from ideal (there are quite a few ways to optimize that particular step), but it's a good starter article to get into curve fitting specifically for Bezier curves. Except for the part I couldn't follow =D – Mike 'Pomax' Kamermans Jun 20 '18 at 1:42