Given a unit vector u in R4, I'd like to build a matrix that performs reflection in the hyperplane orthogonal to u and intersecting the origin.

I considered doing it by conjugating the standard reflection matrix in the hyperplane orthogonal to (1,0,0,0) by a rotation matrix, but I realized I'm not sure how to determine the conjugating rotation matrix. In three dimensions, one could for instance compute the axis to be the cross product of u and (1,0,0) and the angle to be the inverse cosine of the dot product of those two. But rotations in four dimensions aren't so simple and may have an entire 2-plane as an invariant axis or they may consist of two simultaneously rotations in orthogonal 2-planes and hence fix only a single point, something impossible in three dimensions.

Any ideas?

You know how to project onto the plane; you just subtract off the component normal to it. To reflect, you instead subtract off TWICE the component normal to it. That's all.

See Householder reflections.

Thanks very much.