Calculating the determinant of a matrix using a purely analytical method that involves the "cross product" in nD

Question

I've been conversing with the various AIs about the cross product in 3D and 7D. I've been lead to believe that n = 3 and n = 7 are only dimensions where the cross product is a well-defined operator. So I took the 7D code and made it into a template function that calculates for any n > 0, not just n = 3 and n = 7. I was not expecting it to work, but it does. Here I apply it as a method of obtaining the determinant of a square matrix.

My chat log with Claude is at: https://claude.ai/chat/3caf4077-28b5-497f-b704-1b0c336a104d

I asked Claude if this was the Hodge star operator, but Claude says that it is not. Please help: what is the name of this operator?

#include <Eigen/Dense>
using namespace Eigen;
#include <iostream>
#include <vector>
using namespace std;
template<size_t N>
class Vector_nD
{
public:
 std::array<double, N> components;
 // Helper function to get the sign of permutation
 static int permutationSign(const std::array<int, (N - 1)>& perm)
 {
 int inversions = 0;
 for (int i = 0; i < (N - 2); i++)
 for (int j = i + 1; j < (N - 1); j++)
 if (perm[i] > perm[j])
 inversions++;
 return (inversions % 2 == 0) ? 1 : -1;
 }
 Vector_nD(const std::array<double, N>& comps) : components(comps) 
 {
 }
 Vector_nD(void)
 {
 for (size_t i = 0; i < N; i++)
 components[i] = 0;
 }
 double operator[](size_t index) const {
 if (index >= N) throw std::out_of_range("Index out of bounds");
 return components[index];
 }
 // Hodge star operator?
 static Vector_nD cross_product(const std::vector<Vector_nD<N>>& vectors)
 {
 if (vectors.size() != (N - 1))
 throw std::invalid_argument("nD cross product requires exactly (n - 1) vectors");
 std::array<double, N> result;// = { 0, 0, 0, 0, 0, 0, 0 };
 for (size_t i = 0; i < N; i++)
 result[i] = 0.0;
 // These are the indices we'll use for each component calculation
 std::array<int, (N - 1)> baseIndices;// = { 0, 1, 2, 3, 4, 5 };
 for (int i = 0; i < (N - 1); i++)
 baseIndices[i] = i;
 // For each component of the result vector
 for (int k = 0; k < N; k++)
 {
 // Skip k in our calculations - this is equivalent to removing the k-th column
 // For each permutation of the remaining (N - 1) indices
 do
 {
 // Calculate sign of this term
 const int sign = permutationSign(baseIndices);
 // Calculate the product for this permutation
 double product = 1.0;
 for (int i = 0; i < (N - 1); i++)
 {
 const int col = baseIndices[i];
 // Adjust column index if it's past k
 const int actualCol = (col >= k) ? col + 1 : col;
 product *= vectors[i][actualCol];
 }
 result[k] += sign * product;
 } while (std::next_permutation(baseIndices.begin(), baseIndices.end()));
 // Reset indices for next component
 for (int i = 0; i < (N - 1); i++)
 baseIndices[i] = i;
 }
 return Vector_nD(result);
 }
 static double dot_product(const Vector_nD<N>& a, const Vector_nD<N>& b)
 {
 double dot_prod = 0;
 for (size_t i = 0; i < N; i++)
 dot_prod += a[i] * b[i];
 return dot_prod;
 }
};
template <typename size_t N>
double determinant_nxn(const MatrixXd& m)
{
 if (m.cols() != m.rows())
 {
 cout << "Matrix must be square" << endl;
 return 0;
 }
 // We will use this vector later, in the dot product operation
 Vector_nD<N> a_vector;
 for (size_t i = 0; i < N; i++)
 a_vector.components[i] = m(0, i);
 // We will use these (N - 1) vectors later, in the cross product operation
 std::vector<Vector_nD<N>> input_vectors;
 for (size_t i = 1; i < N; i++)
 {
 Vector_nD<N> b_vector;
 for (size_t j = 0; j < N; j++)
 b_vector.components[j] = m(i, j);
 input_vectors.push_back(b_vector);
 }
 // Compute the cross product using (N - 1) vectors
 Vector_nD<N> result = Vector_nD<N>::cross_product(input_vectors);
 // Flip handedness
 for (size_t i = 0; i < result.components.size(); i++)
 if (i % 2 == 1)
 result.components[i] = -result.components[i];
 double det = Vector_nD<N>::dot_product(a_vector, result);
 // These numbers should match
 cout << det << endl;
 cout << m.determinant() << endl << endl;
 return det;
}
int main(int argc, char** argv)
{
 srand(static_cast<unsigned int>(time(0)));
 const size_t N = 11; // Anything larger than 11 takes eons to solve for
 MatrixXd m(N, N);
 for (size_t i = 0; i < N; i++)
 {
 for (size_t j = 0; j < N; j++)
 {
 m(i, j) = rand() / static_cast<double>(RAND_MAX);
 if (rand() % 2 == 0)
 m(i, j) = -m(i, j);
 }
 }
 determinant_nxn<N>(m);
 return 0;
}

• shawn_halayka: Also, if you like, the whole code (with your suggested changes) can be found at: github.com/sjhalayka/cross_product_nD

Answer 1

I'll answer some questions about the maths without looking at your code. The thing you're looking for is the wedge product. It is the nD generalization of the cross product. The wedge product satisfies the following properties:

It is anticommutative and multilinear.

\begin{align} u\wedge v &= -v\wedge u\iff u\wedge u=0 &&\text{Anticommutative}\\ (au_1+bu_2)\wedge v &= a u_1\wedge v+bu_2\wedge v&&\text{Multilinear} \end{align}

The wedge product produces elements that are higher dimensional. A wedge product between two vectors,

$$u\wedge v,$$

produces an area element. Taking the wedge between an area element and another vector produces a volume element: $$u\wedge v\wedge w.$$

Because the wedge product between two of the same vectors is zero, the amount of basis elements in each dimension is highly constrained. In 3D there are 3 vectors, 3 area elements and 1 volume element. If we call n the number of basis vectors and k the dimension of the specific element (e.g. area has k=2 and volume has k=3), we get that there are $$\pmatrix{n\\k}=\frac{n!}{k!(n-k)!}$$

elements of that dimension. The "dimension" of such an element is more commonly called the grade. Here you can see the number of elements for each n and k up to 8:

Table with number of multivectors in each dimension

Note that k=0 refers to scalars.

The wedge product between two vectors produces a 2-vector (an area element). In 3D something nice happens: the number of vectors and 2-vectors is the same! In geometric algebra, k-vectors and (n-k)-vectors behave very similarly. So we can take these 2-vectors and map them to normal vectors. This is exactly what the cross product does. This mapping is done by the Hodge star operator, which is an operator that maps k-vectors to (n-k) vectors and vice versa.

In 7D this gets more complicated and I don't understand this myself, but by looking at a smart subspace, we can again restrict ourselves to a 7D output space which can then be transformed into vectors. More info is on wikipedia and the links there mentioned.

So to summarize: the wedge product is the generalized cross product. It works for any number of dimensions. But the output generally won't generally be something that can be interpreted as a vector. In 3D and 7D we can interpret the output as a vector. So in those cases we can make have an operator that takes in two vectors and spits out a vector, as well as having all the properties that you want from a cross product.

---

In response to your comments I took a quick look at the code. The Claude chat didn't work for me and without explanation it is hard to understand what's going on without fully diving into it.

On Wikipedia I found a generalization of the cross product that seems to match your description: the multilinear generalization of the cross product. It takes in n-1 vectors and spits out another vector. It also satisfies the properties your normally expect of the cross product (I'm quoting wikipedia here)

1. The output is perpendicular to all the input vectors.
2. The magnitude is equal to the paralletope spanned by the input vectors, i.e. the multidimensional version of a parallelogram.
3. All the vectors together are positively oriented (right hand rule for multiple dimensions)

The wedge product between n-1 vectors is identical to this. I'll show it in 3D here. ex, ey ez are just the basis vectors

\begin{align} u\wedge v=&(u_xe_x+u_ye_y+u_ze_z)\wedge(v_xe_x+v_ye_y+v_ze_z)\\ =&u_xv_xe_x\wedge e_x + \color{BLUE}{u_xv_ye_x\wedge e_y} + \color{GREEN}{u_xv_ze_x\wedge e_z}+\\ &\color{BLUE}{u_yv_xe_y\wedge e_x} + u_yv_ye_y\wedge e_y + \color{RED}{u_yv_ze_y\wedge e_z}\\ &\color{GREEN}{u_zv_xe_z\wedge e_x} + \color{RED}{u_zv_ye_z\wedge e_y} + u_zv_ze_z\wedge e_z\\ =&\color{RED}{(u_yv_z-u_zv_y)e_y\wedge e_z}\\ &\color{GREEN}{(u_zv_x-u_xv_z)e_z\wedge e_x}\\ &\color{BLUE}{(u_xv_y-u_yv_x)e_x\wedge e_y} \end{align} If you take the hodge star of this expression, ey^ez will get mapped to the vector ex etc.

Now imagine you did this for n-1 vectors: $$v^{(1)}\wedge v^{(2)}\wedge\cdots v^{(n-1)}$$ What will each term look like? it will look something like $$v^{(1)}_{i_1} v^{(2)}_{i_2}\cdots v^{(n-1)}_{i_{n-1}}e_{i_1}\wedge e_{i_2}\wedge\cdots e_{i_{n-1}}$$ Most of these terms will be zero: if any the indices match, i.e. $$i_1=i_2,$$ the term will vanish. Similar to the terms colored black in 3D. The nonzero terms will have a positive sign or a negative sign. To know which one it is, you can sort the basis vectors in the wedge product until they have a certain order. Each time you swap two neighbouring basis vectors, you get a minus sign. Does this start to sound familiar?

• shawn_halayka: Thank you for your detailed response. I am indeed calculating the hypervolume of the parallogram using the nD cross product, to calculate the matrix determinant. This product operation takes in $(n - 1)$ $n$-vectors. I don't think it's quite the same as a $k$-blade sequence of operators, but thank you for your input in any case.
• AccidentalTaylorExpansion: If you take the wedge product between $n$ vectors having $n$ dimensions, you get a number that represents the determinant. The output is an n-hypervolume, but from the table you can see that there is only one hypervolume and so it can be mapped to a scalar. This scalar is the determinant. Hope this helps again.
• shawn_halayka: That's the thing: this method that I posted works (it's used to calculate the matrix determinant), and it gives nD vectors as a result. I'm not really asking why it works, but rather what is its name? It is not the wedge product. Although the wedge product is equivalent in 3-D, the wedge product in 4-D is entirely different from this method that I posted.
• shawn_halayka: Wow! Thanks for the further information!
• shawn_halayka: I wrote a small C++ tutorial on how to calculate the determinant: github.com/sjhalayka/papers/blob/main/claude_op/…

Answer 2

Initial observations

Your indentation and style looks good. Common mistakes include inconsistency or spurious indentation (or lack thereof) but I did not see any of these issues.

Coupling sign to parity in permutationSign

int inversions = 0;
for (int i = 0; i < (N - 2); i++)
 for (int j = i + 1; j < (N - 1); j++)
 if (perm[i] > perm[j])
 inversions++;
return (inversions % 2 == 0) ? 1 : -1;

Can be simplified to:

int sign = 1;
for (int i = 0; i < (N - 2); i++)
 for (int j = i + 1; j < (N - 1); j++)
 if (perm[i] > perm[j])
 sign = -sign;
return sign;

Take advantage of library functions

Since you are using std::array you should take advantage of the built-in methods that class provides.

Vector_nD(void)
{
 for (size_t i = 0; i < N; i++)
 components[i] = 0;
}

Can be simplified to:

Vector_nD(void)
{
 components.fill(0);
}

Bounds checking

There is also a library function that indexes the array but throws an exception if out of bounds.

double operator[](size_t index) const {
 if (index >= N) throw std::out_of_range("Index out of bounds");
 return components[index];
}

Simplifies to:

double operator[](size_t index) const {
 return components.at(index);
}

• shawn_halayka: Thank you for your time and insight. I finally understand the parity checking function. It basically counts the number of swaps that would be needed if a sorting function (any sorting function) were in use. There's no short cut, although Lehmer coding was mentioned by several AI.
• shawn_halayka: It's worth mentioning that the permutation function is the Levi-Civita symbol. Now I see.

Answer 3

• // Reset indices for next component loop is not necessary. The documentation says:

  Returns true if such a "next permutation" exists; otherwise transforms the range into the lexicographically first permutation (as if by std::sort) and returns false.

  so the indices are already in the desired state.

• (削除) permutationSign does not seem necessary. next_permutation produces lexicographically next, and it is easy to prove that the signs just alternate. (削除ここまで)

• In any case, permutationSign works too hard. Its time complexity is \$O(n^2)\$. An inversion count (what essentially it does) can be computed in \$O(n\log{n})\$. It is not really important for small \$n\$, but still...

• dot_product is std::inner_product (don't forget to #include <numeric>).

• You may want to go an extra mile, and templatize Vector_nD to work with value types other than double.

• cross_product indeed looks like a special case of Hodge star.

• user555045: "it is easy to prove that the signs just alternate" that's what I thought would happen, but then I tried it and it didn't happen, of course it's possible that I did it wrong..
• vnp: @user555045 I stand corrected. There are ways to compute sign while generating next permutation, but of course they rule out std::next_permutation.
• vnp: @shawn_halayka A special case in \$k = n-1\$ sense.
• vnp: @shawn_halayka Correct
• user555045: @shawn_halayka it's possible, but another option is using a different permutation-enumeration algorithm that does have the property that the sign alternates (then you don't need to store anything either). std::next_permutation doesn't have that property as we found out, but Steinhaus–Johnson–Trotter does

Answer 4

This is not a cross product, which refers to an operation on two vectors yielding a vector (binary operation). The 7D special case is a red herring because you are not using the 7D cross product when n = 7. Instead, your operation maps (n - 1) vectors to a vector, which is well-defined in any dimension: It is the Hodge dual of the wedge product.

In physical terms, the wedge product of (n - 1) vectors defines the size and orientation of a hypersurface element, which for your problem is the "base" of a parallelepiped. The Hodge star operator converts this into a (pseudo)vector that is normal (perpendicular) to the hypersurface. Then, your dot product with the remaining vector computes the hypervolume (determinant) of the parallelepiped.

This determinant can also be interpreted directly as the Hodge dual of the wedge product of all n vectors. It is a (pseudo)scalar -- an "oriented" hypervolume that changes sign upon mirror reflection.

Only in 3D does the surface normal vector correspond to a cross product.

• shawn_halayka: Thank you for the insight. So the operator is the Hodge star operator, where k = n - 1?
• nanoman: @shawn_halayka It is the Hodge star composed with the (n - 1)-fold wedge product.
• shawn_halayka: That's awesome to know. Thank you for confirming that it's the Hodge star operator.
• shawn_halayka: So is this the proper way to write the operator? *(a ^ b ^ c ... etc ) ?