I would've derived it differnetly. Rather than explicitly setting up the equation for squared distances then finding the solution that minimizes it, I would've just used the normal equations to do the least squares part for me. IMO this makes it a bit more "mechanical" in figuring out the soluion (I just need to figure out what equation I'm trying to solve, the normal equation can be memorized and will figure out how to compute the least squares solution to it).
So we already have the equation from the stackexchange thing for computing the "offset vector" from a point to the line. You can express this in different ways, but they way they did it is fine. Basically the 3D vector telling us how far a point x is from the line described by an origin o and a direction u.
ray_to_point_offset = (I - u*ut)*x - (I-u*ut)*o
That first term is "x projected onto a plane with normal u", and the second term is "ray origin projected to the same plane". So the difference between the two is just finding the perpendicular offset vector between the ray and x, which is what we're trying to minimize.
So we basically want to set this equation to zero and try to solve for x. In general we won't be able to find a x where this offset vector is exactly zero but it's what we're striving for (and least squares will minimize the squared residual of this equation, which in this case is just the squared distance from x to the ray).
(I - u*ut)*x - (I-u*ut)*o = 0
(I - u*ut)*x = (I-u*ut)*o
A*x = b, where A = (I - u*ut) and b = -(I-u*ut)*o
A is now a 3x3 matrix, and b is a 3x1 matrix (i.e. a vector). In other words, there are three rows in this equation system.
This is just for one ray, to solve for multiple of them, we vertically "stack" these equations for each ray. E.g. for 4 rays, we would have 4 of these 3-row equations above, so we stack vertically to construct a 12 row system. In other words, in the A*x=b form, A would be 12x3 and b would be 12x1. x is still just a 3 element (column) vector representing the point we're trying to find. You can stack as many of these as we want (e.g. you could have a 300x3 "A" matrix if you have 100 rays).
Now this is over-determined (3 unknows, 12 equations), so to get the least squares solution we construct the normal equation (see wikipedia, or just use the formula on trust): At*A*x = At*b
Solving for x in the system above gives you the least squares solution.
At*A is a 3x3 matrix (for 4 rays, A is 12x3, so At is 3x12. Multyply the two and you get 3x3). You can construct this matrix directly rather than first computing a 12x3 matrix... Similary At*b is 3x1 and you can construct that directly too rather than first constructing the 12x1 b vector and multypling it by the 3x12 At matrix.
Solve for x, in this case by inverting (but really it's just a 3 variable equation system, any simple method of solving this will do here):
x = inverse(At*A)*At*b
This amounts to inverting a 3x3 matrix (At*A) and multpyling it by a 3x1 vector (At*b). Note, it's always 3x3 regardless of how many rays you have.
Simple code example:
Mtx3x3 AtA = Mtx3x3::zero();
Vec3 Atb = Vec3::zero();
for (int r = 0; r < _countof(rays); ++r) {
Vec3 dir = rays[r].dir;
Vec3 origin = rays[r].origin;
// Compute the per-ray "A" matrix.. I.e. the projection matrix
Mtx3x3 A = Mtx3x3::identity() - Mtx3x3{
Vec3{ dir.x*dir.x, dir.x*dir.y, dir.x*dir.z},
Vec3{ dir.y*dir.x, dir.y*dir.y, dir.y*dir.z },
Vec3{ dir.z*dir.x, dir.z*dir.y, dir.z*dir.z },
};
// Not really necessary for this particular problem, since A is symmetric.. but for clarity...
Mtx3x3 At = A.transpose();
// Add up our part of AtA for the full system
AtA += At * A;
// Now add up At*b
Vec3 b = A*origin;
Atb += At*b;
}
Vec3 x = AtA.invert()*Atb;