Orthogonalize[ ] working as expected only when applied twice

Question

Applying Orthogonalize[] once:

v1 = PolyhedronData["Dodecahedron", "VertexCoordinates"][[1]]; 
Graphics3D[Line[{{0, 0, 0}, #}] & /@
  Orthogonalize[{a, b, c} /.
    FindInstance[{a, b, c}.v1 == 0 && (Chop@a != 0.||Chop@b != 0.||Chop@c != 0.),
     {a, b, c}, Reals, 4]], Boxed -> False]

And now twice:

Graphics3D[Line[{{0, 0, 0}, #}] & /@
  Orthogonalize@Orthogonalize[{a, b, c} /.
    FindInstance[{a, b, c}.v1 == 0 && (Chop@a != 0.||Chop@b != 0.||Chop@c != 0.),
     {a, b, c}, Reals, 4]], Boxed -> False]

Errr ... Why?

Answer 1

I think the first result is due to numerical error, taking

sys = {a,b,c}/.FindInstance[
          {a, b, c}.v1 == 0 && (Chop@a != 0. || Chop@b != 0. || Chop@c !=0.), 
          {a, b, c}, Reals, 4];

then MatrixRank@sys returns 2, therefor the system itself is only two dimensional. To me, this implies that the first instance of Orthogonalize is generating a numerical error, and the second instance is using the out of plane error to give you your three vectors. Removing the Chop conditions fixes this,

Orthogonalize[{a, b, c} /.
    N@FindInstance[{a, b, c}.v1 == 0,{a, b, c}, Reals, 4]]

where N is necessary to get rid of the Root terms that appear. This gives you a two-dimensional system, but you can get a third by taking the cross product.

Edit: Here's further evidence that its numerical error due to Chop.

With Chop, FindInstance gives me

{{64., 3.6, 335.108}, {-67., -4.3, -350.817}, {0, 176., 0}, 
 {-2., -4.3, -10.4721}}

Without Chop, I get

{{-16.8, 3.9, -87.9659}, {6.6, -1.7, 34.558}, {13.4, -4.3, 70.1633}, 
 {19.9, -4.3, 104.198}}

which is a significant difference between the two.

Answer 2

I also assumed it would be a numerical error, but didn't quite understand why, so I tried to implement Gram-Schmidt orthogonalization myself, hoping to understand the problem on the way:

(* projects onto a unit vector *)
proj[u_][v_] := (u.v) u

Clear[gm, gramSchmidt]

gm[finished_, {next_, rest___}] := 
 With[{v = next - Plus @@ Through[(proj /@ finished)[next]]}, 
  gm[Append[finished, Normalize@Chop[v]], {rest}]
 ]

gm[finished_, {}] := finished

gramSchmidt[vectors_] := gm[{}, vectors]

^{(Included for illustration only, I simply couldn't quite figure out what's going on before I reimplemented it myself.)}

A critical step here, which I didn't realize before, is deciding whether a vector we get is zero or not before the normalization step (see Chop in my code). Otherwise we might get something tiny, possibly a mere numerical error, which is then normalized back into a large value.

This seems to be controlled by the Tolerance option of Orthogonalize, and indeed, raising the tolerance, and forcing it to discard tiny vectors fixes the problem you describe. Orthogonalize[ ... , Tolerance -> 1*^-10] works in a single step.

Answer 3

Perhaps it is a characteristic of the default GramSchmidt method?

Try: Method -> "Reorthogonalization" or Method -> "Householder".