∇ΘJ=∇Θ(y−XΘ)T(y−XΘ)
∇ΘJ=∇ΘyTy−(XΘ)Ty−yTXΘ+ΘT(XTX)Θ
https://news.ycombinator.com/item?id=21301117
in this thread. I keep the math and notation really simple. The only prerequisite is matrix multiplication. I use no calculus or probability. I derive the normal equations in just two simple lines of matrix algebra.
If you want software, there is lots of it or you can write your own easily enough: For the matrix work, might use LINPACK or just write a little Gauss elimination routine for yourself.
It's not, there's an implied parenthesis:
(∇_Θ)J = (∇_Θ)( yTy − (XΘ)Ty − yTXΘ + ΘT(XTX)Θ ).
In the subsequent expressions you can see derivatives being taken of all the terms.
1. Expansion. Let J(x) = norm(Ax-b)^2. Then,
J(x) = norm(Ax-b)^2 = (Ax-b)'(Ax-b) = (x'A'-b')(Ax-b) = x'A'Ax - b'Ax -x'A'b + b'b = x'A'Ax - 2b'Ax + b'b
Then,
grad J(x) = 2A'Ax-2A'b
Except, how did we know that grad(x'A'Ax) = 2A'Ax or grad(2b'Ax)=2A'b? Well, there are at least two ways to derive it. One, go back to the limit definition. This is tedious. Two, cheat and use a Taylor series. Now, this assumes that we believe Taylor's theorem and it's worthwhile to go through a proof someday. Today's not that day. Taylor series assert that under the proper conditions that:
J(x+dx) = J(x) + J'(x)dx + (1/2)J''(x)dxdx + R(x)
where J'(x) denotes the total derivative of J and R(x) is a remainder term. This is a linear operator where if J:X->Y, then J'(x)\in L(X,Y) where L(X,Y) denotes the space of linear operators that have domain X and codomain Y. This means that J''(x)\in L(X,L(X,Y)) since we can still play this same trick. This means that J''(x)dxdx \in Y. Anyway, why do we care? If we can write something that looks like a Taylor series, then we can read off the derivatives. For example, for K(x) = x'Bx, we have that:
K(x+dx) = (x+dx)'B(x+dx) = x'Bx+dx'Bx+x'Bdx+dx'Bdx = x'Bx+2x'Bdx+dx'Bdx = K(x)+K'(x)dx+1/2K''(x)dxdx
Hence, if we just match the terms:
K(x)=x'Bx K'(x)dx=2x'Bdx (1/2)K''(x)dx=dx'Bdx or K''(x)=2dx'Bdx
Great, except grad K(x)=2B'x and not 2x'B. Why the transpose? Well, it's because K'(x) and grad K(x) are not the same thing. Recall, from above, if K:X->R where R means the real numbers, then K'(x)\in L(X,R). To be clear, this is a function that maps X to R. Say X was something like an m dimensional vector space, R(m). Then, K'(x) : R(m)->R. If we want to represent this as a matrix, the matrix would have size R(1,m). This is a row vector. Normally, the gradient is a column vector. Why? Well, most of the time, people view the gradient as a matrix of partial derivatives and they order them in a magical way. To me, a better way to visualize it is as what's called the Riesz representative of the total derivative. Basically, there's a theorem that states that under certain conditions any linear functional can be represented as the inner product between an element in that space and the argument. Basically, if f is a linear function that maps to the real numbers, then we can write any linear function, under the right conditions, as f(x)=<a,x> where <.,.> denotes inner products. In R(m), this inner product is just the transpose, so any linear functional in R(m) can be written as f(x)=a'x. Going back to the above, K'(x)dx=2x'Bdx=(2B'x)'dx. This means that the linear functional C(dx)=K'(x)dx can be written as C(dx)=K'(x)dx=(2B'x)'dx. That element, 2B'x is the gradient, by definition. It's the Riesz representative of the total derivative. And, again, you can define the gradient differently. I wouldn't since this is clean and generalizes to Hilbert spaces, which are complete inner products spaces, which means we get up to some restriction of infinite dimensions and function spaces.
OK, that was a lot, but then we can realize that we can just apply the Taylor series trick to the whole expansion:
J(x+dx) = norm(A(x+dx)-b)^2 = norm(Ax+b+Adx)^2 = norm(Ax+b)^2 + 2(Ax+b)'Adx + dx'A'Adx = J(x) + J'(x)dx + (1/2)J''(x)dxdx
After matching, we have that
J'(x)dx = 2(Ax+b)'Adx
or that
grad J(x) = 2A'(b+Ax) = 2A'b + 2A'Ax
2. But, wait, there's more. We could just use the chain rule. Recall, (f o g)'(x) = f'(g(x))g'(x) where o denotes function composition. In the above case, f(x)=norm(x)^2 and g(x)=Ax+b. Also playing Taylor series games, we'd eventually find out that f'(x)dx=2 dx' where I is the identity and g'(x)=A. This means that
f'(g(x))g'(x) = 2(Ax-b)'A
which means the gradient is:
2A'(Ax-b)
3. If you really don't like this, you could pretend to be a physicist and do everything component wise with 1-D derivatives. That was a joke. This will produce the correct result, but one that I find particularly painful to derive. As such, I will punt and avoid that here.
Anyway, all teasing aside, there are lots of different ways to get the same result. Really, just use what works well for you and potentially adapt if you're working with an audience that views things differently. I don't claim what I wrote above is the best notation or derivations for you, but they work well for me and I wanted to provide some variety on a Saturday afternoon.