This post is a further development of my earlier question in reply to John's post. I have implemented a basic version of the CANDECOMP/PARAFAC algorithm referred to on Wikipedia and described here. The implementation is not particularly quick or elegant, but if you'd like to experiment with it, you're welcome to copy and paste the code from the startup region of the attached worksheet. (Especially the stopping criterion could use work.)

First let me show you some results. Of course I used my own photo, not John's.

First of all, what do those captions mean? SVD stands for John's method, concatenating the Y, U, and V versions of the image along one of the axes, taking the SVD, then recomposing the result; C/P stands for the CanDecomp/PARAFAC algorithm, viewing the image as a third order tensor and working directly with that. What the rank means should be clear to who has read John's post.

Then it says either "singly-centered" or "multi-centered". I added a little bit of functionality to centre and scale the different "slices" of the third order tensor: centering means subtracting the average, and scaling means making the variance equal to one. This was suggested by the CanDecomp tutorial linked to in the introduction above, and indeed it improved the results. The tutorial claims you should do one such centering for each of the one-dimensional "slices" in one direction through your tensor, because the algorithm doesn't like it if you do it otherwise. So I did that for the CanDecomp algorithm (taking the minimum number of 3*256=768 averages). I did it for the SVD as well, but there only with the three whole colour layers. Scaling should happen with full two-dimensional "slices" at once, so in both cases I scaled each of the three colour slices (Y, U, V) to have variance 1. Then I performed the decomposition and undid first the scaling, then the centering.

With a given rank, the C/P algorithm leads to a stronger "compression": fewer scalars are necessary. This can be seen by considering how a rank-one tensor is specified. For C/P, we have a "singular" vector for the height (with 304 elements in this case), one for the width (256 scalars), and the colours (3 scalars). Each of these is constrained to having norm 1, so the number of degrees of freedom is one less. Adding the singular value itself, we get a total of 303+255+2+1=561 scalars per rank. For the SVD approach, you need a vector for the width (256 scalars) and a vector that consists of three "copies" for the height, one for each colour channel (3*304=912). Again, these two vectors are constrained to have norm one, and we need one more degree of freedom for the singular value, so we need 255+911+1=1167 degrees of freedom per rank for the SVD. To these numbers, we need to add the information for centering (either 3 scalars or 256*3 scalars) and scaling (3 scalars), to obtain the numbers you see above. So you see that the rank 7 C/P image is compressed about as much as the rank 4 SVD image, and I would say it has slightly better image quality (I especially like the eyes).

Finally, of course the C/P algorithm is much slower than SVD, taking as long as 10 seconds to compute the miserable rank 4 approximation you see above on my Core i7 laptop. This is not just because the algorithm is intrinsically more complicated (which it is): my implementation is also just not as efficient. I did try to not allocate new memory for rtables too often inside the inner loops, preallocating memory beforehand instead, most of the time. However, this was not based on any profiling (which I haven't done); I'm sure it could be sped up.

Finally, here is the (very spartan) worksheet I used to obtain the images above: candecomp.mw. That is, this is what I used for the C/P images; the SVD ones are from John's code.

(Edit: removed a typo from the worksheet.)