## Color: Cohen’s intriguing F matrix again

Minor edits, 4 Jan 2010.

Let me run thru the derivation of Figure 20 again, still with the CIE 1964 tables, but at 20 nm intervals. Then I show two alternative calculations. And one of them will show us why Cohen switched the sign on f2.

There is interesting linear algebra and vector algebra in this post, even if you’re not interested in where these problems came from.

## Review:

What we are aiming for is an orthonormal basis whose three vectors are the columns of an F matrix. In particular, the first column — the first basis vector — is the fundamental of an equal-energy illuminant. The challenge was to find two other columns of F which matched Cohen’s Figure 20.

There were scaling issues all over the place. In particular, Cohen’s final three vectors were orthogonal but not orthonormal, because they were not of unit length. Having matched Figure 20, I will no longer concern myself with his scaling; unit vectors, I want unit vectors.

Let’s grab things as we need them.

Since I’ve shown you pictures derived from the 1964 xyz bar tables, I’ll stay with 1964.

What range of wavelengths? The CIE 1964 xyz bar tables run from 380 to 780. What the heck, just use 400 to 700 as before.

What interval? 20 nm, so I can display things.

Then my A matrix is the following subset of the CIE 1964 xyz bar tables:

$A = \left(\begin{array}{ccc} 0.0191 & 0.002 & 0.086 \\ 0.2045 & 0.0214 & 0.9725 \\ 0.3837 & 0.0621 & 1.9673 \\ 0.3023 & 0.1282 & 1.7454 \\ 0.0805 & 0.2536 & 0.7721 \\ 0.0038 & 0.4608 & 0.2185 \\ 0.1177 & 0.7618 & 0.0607 \\ 0.3768 & 0.962 & 0.0137 \\ 0.7052 & 0.9973 & 0. \\ 1.0142 & 0.8689 & 0. \\ 1.124 & 0.6583 & 0. \\ 0.8563 & 0.3981 & 0. \\ 0.4316 & 0.1798 & 0. \\ 0.1526 & 0.0603 & 0. \\ 0.0409 & 0.0159 & 0. \\ 0.0096 & 0.0037 & 0.\end{array}\right)$

There are 16 distinct wavelengths. I take the equal-energy illuminant to have the value 100 for 16 values:

ee= {100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100}

To get its fundamental, I need a projection operator onto the fundamental color space (the pre-image of the range of A’, the transpose of A). To get that, I need an orthonormal basis for the fundamental color space. To get that, I need the Singular Value Decomposition (SVD) of A’:

A’ = u w v’

We know that the first 3 columns (because A only had 3 columns) of v are an orthonormal basis for the fundamental color space; the remaining 13 (= 16 – 3) columns of v are an orthonormal basis for the null space of A’. Here are the first 3 columns of v, and I name them v1:

$v1 = \left(\begin{array}{ccc} -0.0253683 & 0.0153138 & 0.00589786 \\ -0.285062 & 0.175889 & 0.0585877 \\ -0.574584 & 0.357498 & 0.0807567 \\ -0.513035 & 0.307203 & -0.0134903 \\ -0.242873 & 0.0986244 & -0.196418 \\ -0.112652 & -0.0638886 & -0.364177 \\ -0.124169 & -0.198387 & -0.487962 \\ -0.173395 & -0.313102 & -0.431359 \\ -0.221517 & -0.394131 & -0.198767 \\ -0.250764 & -0.427605 & 0.140947 \\ -0.241369 & -0.398876 & 0.386199 \\ -0.171479 & -0.278411 & 0.372286 \\ -0.0839289 & -0.135192 & 0.203385 \\ -0.0292822 & -0.046994 & 0.0743802 \\ -0.00781685 & -0.012531 & 0.020133 \\ -0.00183092 & -0.00293336 & 0.00474977\end{array}\right)$

We know that we construct a projection operator R onto the column space of v1 by

R = v1 v1′.

We can quickly check that R is, at least, a projection operator, since it is idempotent: i.e. R^2 = R.

And we can check that it is a projection operator onto the column space of v1, since it preserves each of the columns of v1: i.e. R v1 = v1.

If we wanted, we could confirm that it projects all of the other columns of v to the zero vector: i.e. R v2 = 0, where v2 is the last 13 columns of v.

Enough already, we came here to get the fundamental of the equal energy illuminant…

f = R ee = {5.53946, 62.01, 125.887, 116.934, 68.0773, 55.4377, 80.9502, 109.171, 126.547, 128.189, 113.072, 76.3005, 36.473, 12.5844, 3.34797, 0.782783}

Make it a unit vector f1…

f1 = {0.0165427, 0.185182, 0.375939, 0.349205, 0.203302, 0.165555, 0.241744, 0.326021, 0.377912, 0.382816, 0.337672, 0.227859, 0.108921, 0.0375811, 0.00999815, 0.00233765}

For the other two vectors, we found that we could match Cohen by taking the three vectors {f1, E2, E3} and applying the Gram-Schmidt orthogonalization procedure… which is nothing more than subtracting off the part of E2 which is parallel to f1 (and normalizing, to get a unit vector f2), and then subtracting off the part of E3 which is in the span of f1 and f2.

First, we need the E basis; they are the dual basis to A’; or E’ is the dual basis to A. As always, I find the solution by writing the definition…

E’A = I

then writing the false (because A’ (correction, A) is not invertible)…

$E' = A^{-1}\$,

which leads me to write the true

$E' = (A'A)^{-1} A'\$.

As usual, that formula is for the transpose E’ but I write E = E”, the transpose of the transpose:

$E = \left(\begin{array}{ccc} 0.00507367 & -0.00518039 & 0.00960814 \\ 0.050171 & -0.0533474 & 0.109524 \\ 0.0711671 & -0.0800759 & 0.224744 \\ -0.00118297 & -0.00392724 & 0.205467 \\ -0.140011 & 0.15312 & 0.103269 \\ -0.256404 & 0.304158 & 0.0467504 \\ -0.323667 & 0.432104 & 0.0296721 \\ -0.247916 & 0.423527 & 0.0110232 \\ -0.0411424 & 0.27367 & -0.0156973 \\ 0.237054 & 0.0289735 & -0.0458278 \\ 0.422342 & -0.164367 & -0.0633192 \\ 0.375936 & -0.19191 & -0.0524917 \\ 0.200409 & -0.110177 & -0.027315 \\ 0.0725721 & -0.0410648 & -0.00979225 \\ 0.0195879 & -0.011175 & -0.00263529 \\ 0.00461444 & -0.00264363 & -0.000619869\end{array}\right)$

Just for my sanity, let me make unit vectors e2 and e3 out of E2 and E3 (which are, in case you have any doubt, the 2nd and 3rd columns of E):

e2 = {-0.00645091, -0.0664311, -0.0997151, -0.00489041, 0.190674, 0.378755, 0.53808, 0.5274, 0.34079, 0.0360794, -0.204679, -0.238978, -0.137198, -0.0511362, -0.0139157, -0.003292}

e3 = {0.02679, 0.305381, 0.626647, 0.572896, 0.28794, 0.130352, 0.0827338, 0.0307357, -0.0437682, -0.12778, -0.176551, -0.146361, -0.0761614, -0.0273034, -0.00734788, -0.00172836}

The part of e2 which is parallel to f1 is…

f1 (e2.f1)
= {0.00585458, 0.0655375, 0.133048, 0.123587, 0.0719501, 0.0585914, 0.0855552, 0.115382, 0.133746, 0.135482, 0.119505, 0.080641, 0.0385479, 0.0133002, 0.00353843, 0.000827313}

so the part of e2 which is orthogonal to f1 is found by subtracting the parallel part from e2…

t1 = e2-f1 (e2.f1)
= {-0.0123055, -0.131969, -0.232763, -0.128477, 0.118724, 0.320164, 0.452525, 0.412018, 0.207044, -0.0994022, -0.324183, -0.319619, -0.175746, -0.0644364, -0.0174541, -0.00411931}

Do a quick check that t1 is othogonal to f1…

t1.f1 = -5.15775*10^-17

But t1 is not a unit vector: it’s squared length is t1.t1 = 0.874749 .

So we normalize it, and call the result f2:

f2 = {-0.013157, -0.141101, -0.24887, -0.137367, 0.126939, 0.342318, 0.483839, 0.440529, 0.221371, -0.106281, -0.346616, -0.341736, -0.187908, -0.0688953, -0.0186619, -0.00440436}

To get f3, we find the part of e3 which lies in the plane spanned by f1 and f2…

f1 (f1.e3)+f2 (f2.e3)
= {0.00735195, 0.0822258, 0.166479, 0.153523, 0.0869147, 0.0679319, 0.0993851, 0.136557, 0.161737, 0.167772, 0.150999, 0.103177, 0.0496131, 0.0171664, 0.00457096, 0.00106922}

then subtract that from e3, and the result will be orthogonal to the plane spanned by f1 and f2…

t2 = e3-(f1 (f1.e3)+f2 (f2.e3))
= {0.0194381, 0.223156, 0.460167, 0.419372, 0.201026, 0.0624205, -0.0166512, -0.105822, -0.205505, -0.295552, -0.32755, -0.249538, -0.125774, -0.0444698, -0.0119188, -0.00279758}

Like t1, t2 is not a unit vector, but we make it so, and call the result f3:

f3 = {0.0215889, 0.247847, 0.511084, 0.465775, 0.223269, 0.0693272, -0.0184937, -0.117531, -0.228244, -0.328255, -0.363793, -0.277149, -0.139691, -0.0493903, -0.0132376, -0.00310713}

If we arrange f1, f2, f3 as the columns of F…

$F = \left(\begin{array}{ccc} 0.0165427 & -0.013157 & 0.0215889 \\ 0.185182 & -0.141101 & 0.247847 \\ 0.375939 & -0.24887 & 0.511084 \\ 0.349205 & -0.137367 & 0.465775 \\ 0.203302 & 0.126939 & 0.223269 \\ 0.165555 & 0.342318 & 0.0693272 \\ 0.241744 & 0.483839 & -0.0184937 \\ 0.326021 & 0.440529 & -0.117531 \\ 0.377912 & 0.221371 & -0.228244 \\ 0.382816 & -0.106281 & -0.328255 \\ 0.337672 & -0.346616 & -0.363793 \\ 0.227859 & -0.341736 & -0.277149 \\ 0.108921 & -0.187908 & -0.139691 \\ 0.0375811 & -0.0688953 & -0.0493903 \\ 0.00999815 & -0.0186619 & -0.0132376 \\ 0.00233765 & -0.00440436 & -0.00310713\end{array}\right)$

we can compute F’F. Since F is supposed to be orthonormal, we expect to get a 3×3 identity matrix. Sure enough…

$F'F = \left(\begin{array}{ccc} 1. & 0 & 0 \\ 0 & 1. & 0 \\ 0 & 0 & 1.\end{array}\right)$

The resulting F, of course, are far from smooth…

but they still have the same kind of shapes we saw in Figure 20.

Oh, except that to match his figure, I had to take the negative of f2.

But I know of no good reason to do that, so I’m not going to. Yet. (Ok, Ok — the obvious guess is that f1, f2, f3 is a left-handed basis. If we change the sign of f2, we will get a right-handed basis. I’ll show it to you later in this post. That is, there is a good reason to do that, but nothing we’ve seen yet would lead us to it. )

What have we done? On a coarser scale, we have reproduced the calculations that lead to a specific orthonormal F matrix. It is one of many — infinitely many — bases which have f1 as the first basis vector. That really is the only special vector; the other two span the plane perpendicular to f1.

## Cohen’s Cholesky Decomposition

Now let’s try it the way I think Cohen did it.

There is one potentially confusing problem: notation, as so often happens.

In the following discussion, “A” is no longer a matrix of color-matching functions. We will take “A” to have columns f, E2, E3; that was my starting point for the Gram-Schmidt orthogonalization. He describes such a starting point, but I believe he never shows the computations for such a case.

Here then is our new A matrix (which I call A1 in practice):

$A1 = \left(\begin{array}{ccc} 5.53946 & -0.00518039 & 0.00960814 \\ 62.01 & -0.0533474 & 0.109524 \\ 125.887 & -0.0800759 & 0.224744 \\ 116.934 & -0.00392724 & 0.205467 \\ 68.0773 & 0.15312 & 0.103269 \\ 55.4377 & 0.304158 & 0.0467504 \\ 80.9502 & 0.432104 & 0.0296721 \\ 109.171 & 0.423527 & 0.0110232 \\ 126.547 & 0.27367 & -0.0156973 \\ 128.189 & 0.0289735 & -0.0458278 \\ 113.072 & -0.164367 & -0.0633192 \\ 76.3005 & -0.19191 & -0.0524917 \\ 36.473 & -0.110177 & -0.027315 \\ 12.5844 & -0.0410648 & -0.00979225 \\ 3.34797 & -0.011175 & -0.00263529 \\ 0.782783 & -0.00264363 & -0.000619869\end{array}\right)$

Cohen computed two F matrices (F1 and F2) using two Cholesky decompositions, but we only need one of them. We recall that the first column of F2 = E Ga is proportional to the first column of A — and that’s what we want — so we will need E1, i.e. the dual basis to A1.

That is crucial. We need the dual basis to A1, not to A. We’ve already used the dual basis to A, to get two columns of A1, but now we need the analog of E for A1 instead of for A.

By definition, (dropping the 1’s for a moment)

E’A = I

and so

$E' = (A'A)^{-1} A' = Me\ A'\$.

I am using Cohen’s

Ma = A’A

Me = E’E

and it turns out that Ma and Me are inverses. I will, of course, apply that equation for A and E to A1 and E1.

Here is Ma and its inverse Me:

$Ma = \left(\begin{array}{ccc} 112130. & 95.1685 & 52.2359 \\ 95.1685 & 0.644885 & 0.0411263 \\ 52.2359 & 0.0411263 & 0.128627\end{array}\right)$

$Me = \left(\begin{array}{ccc} 0.0000123197 & -0.00153021 & -0.00451381 \\ -0.00153021 & 1.77301 & 0.054536 \\ -0.00451381 & 0.054536 & 9.59006\end{array}\right)$

Here is E1:

Now we do a Cholesky decomposition of Ma:

$Ga = \left(\begin{array}{ccc} 334.859 & 0. & 0. \\ 0.284205 & 0.751074 & 0. \\ 0.155994 & -0.00427115 & 0.322916\end{array}\right)$

That is, the Cholesky decomposition has given us Ga such that Ga Ga’ = Ma.

Cohen had showed us that F2 = E Ga was an orthonormal matrix:

$F2 = E Ga = \left(\begin{array}{ccc} 0.0165427 & -0.013157 & 0.0215889 \\ 0.185182 & -0.141101 & 0.247847 \\ 0.375939 & -0.24887 & 0.511084 \\ 0.349205 & -0.137367 & 0.465775 \\ 0.203302 & 0.126939 & 0.223269 \\ 0.165555 & 0.342318 & 0.0693272 \\ 0.241744 & 0.483839 & -0.0184937 \\ 0.326021 & 0.440529 & -0.117531 \\ 0.377912 & 0.221371 & -0.228244 \\ 0.382816 & -0.106281 & -0.328255 \\ 0.337672 & -0.346616 & -0.363793 \\ 0.227859 & -0.341736 & -0.277149 \\ 0.108921 & -0.187908 & -0.139691 \\ 0.0375811 & -0.0688953 & -0.0493903 \\ 0.00999815 & -0.0186619 & -0.0132376 \\ 0.00233765 & -0.00440436 & -0.00310713\end{array}\right)$

I want to point out now, that the lower-triangular structure of Ga says that the 3rd column of F2 is proportional to the third column of E1: in fact, I claim we have (although it’s difficult to make it clear)

$F2_3 = .322916 E1_3$

(i.e. the 3rd column of F2 is .322916 times the 3rd column of E1).

Cohen also told us (I’ll confess that I haven’t figured out why) that the 1st column of F2 is proportional to the first column of A1: in fact, we have

$F2_1 = .00298633 A1_1$

Anyway, we have just applied his method to (f, E2, E3). And we have confirmed the realtionships between F2 and A1 and E1.

Fine, but what did we get? Did we get the same answers as from Gram-Schmidt?

Yes. The first column of F2 is equal to f1, the second column to f2, and the third column of F2 is equal to f3. (I could put up a drawing… but it looks exactly like the previous one.)

Lovely! And I’ll point out that we matched my f2 — not the negative of it. Even his method would require that we take the negative of his answer in order to get Figure 20.

## The cross product

There is one other alternative computation I want to do.

We just saw that f1 was proportional to the first column of A1 — that’s how we chose it, a unit vector proportional to the fundamental f.

But we also saw that f3 was proportional to the third column of E1 — and so we could have computed f3 as soon as we had E1, just by normalizing the third column of E1. Cohen seems to have gone to an awful lot of work to get the matrix F2, given that both its first and third columns are known as soon as we compute E1 from A1. (This would work in all of his cases: for F2, normalize the first column of A and the third column of E; or — for F1 — normalize the first column of E and the third column of A. Edit: In either case, take the first column of one and the third column of the other.)

Anyway, that leaves f2, the second column of F2, to be determined.

Gee, it’s too bad we can’t compute a cross product for vectors of length 16.

Gee, who needs to?

Those vectors lie in a 3D subspace, so there is a change-of-basis transition matrix which leaves them with 3 nonzero components. Edit: and in 3D, we can indeed take the cross product of two vectors to get a third one orthogonal to them. f2 can be computed from f1 and f3 by a cross product. In fact, the transition matrix is precisely v, and we know that new components y are related to old components x by

x = v y

so

$y = v^{-1} x = v'x\$

since v is orthogonal. For old components f1, we get

t1 = v’ f1
= {-0.913557, -0.393456, -0.102984, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}

so we take y1 to be the three nonzero components…

y1 = {-0.913557, -0.393456, -0.102984}

For old components f3, we get

t3 = v’ f3
= {-0.361692, 0.901754, -0.236681, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}

and, again, we take y3 to be the three nonzero components…

y3 = {-0.361692, 0.901754, -0.236681}

Now, what we have are the i and k unit vectors (first and third columns), we want the j unit vector (second column). The right-hand-rule shows me — you can’t do this with your right hand on the mouse! — that

j = k x i,

so I compute

$y2 = y3 \times y1$

y2 = {-0.18599, 0.178973, 0.966114}

(which is a unit vector) and use

x = v y

to get the corresponding f. Oh, I need to pad y2 out with zeroes.

t4 = {-0.18599, 0.178973, 0.966114, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}

t2= v t4
= {0.013157, 0.141101, 0.24887, 0.137367, -0.126939, -0.342318, -0.483839, -0.440529, -0.221371, 0.106281, 0.346616, 0.341736, 0.187908, 0.0688953, 0.0186619, 0.00440436}

Quickly recall f2

f2 = {-0.013157, -0.141101, -0.24887, -0.137367, 0.126939, 0.342318, 0.483839, 0.440529, 0.221371, -0.106281, -0.346616, -0.34,1736, -0.187908, -0.0688953, -0.0186619, -0.00440436}

and we can confirm that t2 = – f2.

and that — presumably — is why he took the negative of f2: both his method and my Gram-Schmidt delivered left-handed coordinate systems; the cross product let us construct a right-handed coordinate system.

I think that’s pretty slick. Sure, it’s just a cross-product, but it was in a 3D subspace of $R^{16}\$.