Cracovians: The Twisted Twins of Matrices

Linear algebra is typically explained using matrices. But matrix theory is just one possible perspective. Below, I describe an alternative approach to linear algebra.

Tadeusz Banachiewicz (1882–1954), a Polish astronomer living in Krakow, was passionate about calculating machines. From the 1920s, Banachiewicz developed a method for computations on tables of numbers, which was particularly easy to perform with arithmometers. In honor of Krakow, Banachiewicz named these computational objects cracovians. In the preface to Banachiewicz’s book, Cracovian Algebra (Rachunek krakowianowy, PWN 1959), it is noted that “some calculations, generally known to be very tedious, became a sheer entertainment for the executor.”

Similar to matrices, cracovians are represented as rectangular tables of numbers. The equality of cracovians, addition of cracovians, and multiplication of cracovians by a scalar are defined identically to their matrix counterparts. However, the multiplication of a cracovian by another cracovian is defined differently: the result of multiplying an element from column i of the left cracovian by an element from column j of the right cracovian is a term of the sum in column i and row j of the result. In essence, each column of the left cracovian is multiplied by each column of the right cracovian.

Here’s an example (cracovians are enclosed in braces to distinguish them from matrices):

All such cracovians whose elements on the main diagonal are equal to 1 and other elements are zeros are denoted by the Greek letter τ (tau). Any such unit cracovian τ as the second factor of multiplication does not change the cracovian by which it is multiplied. As the first factor of multiplication, a unit cracovian τ transposes the cracovian by which it is multiplied:

For any such cracovians A and B that have the same number of rows, the identity holds:

Therefore, in the general case, cracovian multiplication is not commutative.

By convention, we agree that cracovian multiplication is left-associative, i.e.,

The product of cracovians has as many columns as the first factor, and as many rows as the last factor has columns.

The following identities hold:

and so on. All products that occur when calculating τ(A₁​×…×A_n​) have the same number of columns as A_n. ​If the last factor A_n​ has fewer columns than the first factor A₁​, then it is simpler to calculate τ(A₁​×…×A_n) than A₁​×…×A_n.

The following identities hold:

Thus, in the general case, cracovian multiplication is not associative, meaning that cracovians with multiplication do not form a group.

By decomposing a cracovian into a product of two upper triangular cracovians, we can solve systems of linear equations, for example:

which is

In general:

When we decompose the cracovian P into a product of two upper triangular cracovians GH, the following cracovian equations are equivalent:

Since

we can compute z, y, and x:

Banachiewicz rediscovered Cholesky decomposition: any positive-definite symmetric cracovian can be square-rooted, resulting in an upper triangular cracovian. He presented two algorithms for solving symmetric systems of equations.

Furthermore, Banachiewicz presented a cracovian approach to the method of least squares and, among other things, the following applications of quaternions:

In algebra:

• Continued fraction reduction
• Polynomial division
• Horner’s scheme
• Interpolation of functions of two arguments directly from function values, without creating differences
• Iterative solution of systems of linear equations
• Determination of the root of an algebraic equation with the smallest absolute value
• Determination of the greatest common divisor of numbers in a sequence whose terms are given by a polynomial with integer coefficients

In astronomy:

• Determination of the cosines of six angles between the orbital axis directions and coordinate axes from orbital elements or from two heliocentric positions
• Determination of topocentric lunar coordinates
• Differential corrections for planetary and cometary orbits
• Conformal mapping of a spheroid onto a plane
• Determination of precessional coefficients

In geodesy:

• Forward intersection
• Resection
• Hansen’s problem

In spherical trigonometry:

• Composition and decomposition of rotations
• Relationships between elements of a spherical polygon
• Gauss’s formulas for spherical trigonometry
• Solution of spherical hexagons, pentagons, and quadrilaterals
• Quaternion composition and decomposition of rotations

In Banachiewicz’s book, one can find charming advice, such as this: “In practice, when multiplying cracovians, we use movable objects (pen nibs, pencils, matches, etc.) to help the eye in selecting elements from the i-th and j-th columns of the multiplied cracovians. We place these objects on paper at the head (or side) of the appropriate column.”

In 2025, people no longer multiply large matrices by hand. It turns out that multiplying cracovians by computers is not faster than multiplying matrices.

NumPy matrices have two special fields:

• shape: Defines the number of elements in each dimension.
• strides: Determine how many bytes you need to jump in memory to get to the next element in a given dimension.

When NumPy transposes a matrix, it doesn’t create a new matrix with copied data. Instead, NumPy creates a “view” on the same data in memory. This view preserves the contents of the old matrix and swaps the two elements of its shape field and the two elements of its strides field.

Using the timeit.repeat function, we observe that the time for matrix multiplication AB is similar to the time for matrix multiplication B^TA, which is equivalent to the cracovian multiplication AB.

!/usr/bin/env python3

import timeit

import numpy as np

N = 4096

def SetUp():

    global A, B, BT

    A = np.random.rand(N, N).astype(np.float32)

    B = np.random.rand(N, N).astype(np.float32)

    BT = np.transpose(B)

print(timeit.repeat(lambda: np.matmul(A, B), setup=SetUp, number=1))

print(timeit.repeat(lambda: np.matmul(BT, A), setup=SetUp, number=1))

print(A.shape, A.strides, B.shape, B.strides, BT.shape, BT.strides)

Here is the output of the program:

[0.35147937503643334, 0.36795954196713865, 0.3467264170758426, 0.32935020816512406, 0.3452573330141604]

[0.3460142908152193, 0.3827329999767244, 0.3319119580555707, 0.3326827920973301, 0.32132083317264915]

(4096, 4096) (16384, 4) (4096, 4096) (16384, 4) (4096, 4096) (4, 16384)