Finding the most orthogonal set of n vectors from dataset of unit vectors

I am attempting to find the set of determine a subset of vectors in a matrix that are the most dissimilar. Each row in this matrix is a unit vector which is unique.

idx = nchoosek(1:size(dataset, 1), 5);  % Produce all possible combinations of indices
score = zeros(size(idx, 1), 1);  % Store similarity score of unit vectors into the corresponding row of score
for I = 1 : size(idx, 1)
    score(I) = sum(prod(dataset(idx(I, :), :), 1));  % Element-wise multiply this combination's vectors and add
end
[best_scores, best_idx] = mink(score, 10)

[David Goodmanson]

Hi Alex

cosider

a = [ 1 1 1 1 1 1 1 1 -2

2 2 2 2 2 2 2 2 -4

3 3 3 3 3 3 3 3 -6]

Since these three vectors are scalar multiples of each other, they are parallel and as mutually nonorthogonal as you can get. But sum(prod (a)) = 0. So that test is not goiing to work.

[the use of 1 in prod(a,1) is superfluous since 1 is the default]

Answer 1

[Bruno Luong]

I propose this

% random test A, which represents 100 normalized vector in R^5
n = 5;
m = 100;
A = randn(n,m);
A = A ./ sqrt(sum(A.^2,1));

% Want to select k columns of A, most mutually "orthogonal"
k = 3; % requirement: k <= rank(A) <= n
n = size(A,2);
c = nchoosek(1:n,k);

d = arrayfun(@(r) detsub(A(:,c(r,:))), 1:size(c,1));
[~,imax] = max(d);

% selection of best k columns of A
col = c(imax,:);
B = A(:,col)
% check scalar products, it must have small off-diagonal terms
B'*B

function d = detsub(A)
[~,R] = qr(A,0);
d =prod(abs(diag(R)));
end

[Z]

Please let me ask if I am understanding this correctly. Is this effectively taking a subset of the column vectors, forming a set of orthonormal basis vectors from them, then taking the product of all the dot products of the subset vectors with it's accompanying unit vector in the QR decomposition? Does this then mean that the maximal product of dot products for a subset corresponds to that subset's vectors being the best match to it's orthonormalization and therefore the most inherently orthogonal combination of k vectors from the original matrix?

[Z]

If this is correct, then may I ask a question: Why do you judge the subset's match to it's orthogonalization according to the product of it's dot products with the orthonormal basis vectors instead of the sum of dot products of the subset's vectors with the orthonormal basis vectors?

[Bruno Luong]

I don't endorse your interpretation. However you are free to interpret that way if it helps your intuition. But , no, mathematically your interpretation goes far beyond the calculation of determinant suggests.

[Z]

Forgive me, but I don't see where a determinant is being computed in your example code. I see the the QR decomposition is computed on a given subset of column vectors, then the product of the diagonal elements of the R matrix is returned as the metric of orthogonality for that subset. Based on my reading, aren't the diagonals of the R matrix the inner products of the input matrix's column vectors with the orthonormal basis Q computed by the QR decomposition?

[Bruno Luong]

This is the (absolute value) determinant of A (projected on span { A(:,j) } which is span { Q(:,j) )}

prod(abs(diag(R)))

[Z]

Ah, so the absolute value determinant of A(:, j) is can be computed by taking its projection onto it's corresponding orthonormalization Q which represents a unit n-cube, and taking the product of the diagonal elements of R is a way to efficienty compute that. A less ideally and inherently orthogonal A(:, j) set of basis vectors would have a smaller score as measured by this projection onto Q, and a better A(:, j) would have a larger projection / volume according to Q. Is this a more accurate description?

[Bruno Luong]

Yes

Let me explain to you more. When you do qr on A where each column of A is unit vector:

[Q,R] = qr(A)

you have the following properties

• span {Q(:,1:i)} contains span {A(:,1:i)} for all i.
• the coefficient R(i,j) is dot(Q(:,i),A(:,j)) for i <= j.
• Since norm(A(:,j))=1 and R is triangular we have norm(R(1:j,j)) = 1.

Now doing some math. For i < j.

abs(dot(A(:,i),A(:j)) =  norm Projection A(:,j) on span{A(:,i)}
                      <=  norm Projection A(:,j) on span{Q(:,1:i)} 
                      = norm(R(1:i;j))
                      = sqrt(R(1,j)^2 + R(2,j)^2 + ... R(i,j)^2)                      

You want abs(dot(A(:,i),A(:j)) small (most orthogonal) meanning that you want

                  R(1,j)^2 + R(2,j)^2 + ... R(i,j)^2

to be small, for all i < j.

Since from the third property

              R(1,j)^2 + R(2,j)^2 + ... R(j,j)^2 = 1

you just want the last term R(j,j)^2 is as large as possible (therefore the partial sum is small).

And that you want for all j.

I chose det(A) = prod(R(j,j)) as metric, since it has some geometrical interpretation as the volume of the parallelepipied formed by columns of A.

Argually you could use other metric such as "Trace"(A) = sum (abs (R(j,j)) as metric.

My own intuition tell me the det choice is somewhat better.

[Z]

That is a valuable explanation. Thank you so much for your input!