Compression of Multivariate Phase Signals
Your signal is easily described by two parameters (the initial value and the frequency), the sawtooth may be linearized, by stitching each linear segment to the end of the other, and then you can use a linear regression.
The function that stitch the linear segments of a sawtooth of amplitude $2 \pi$ is the unwrap, there is a function with the same name in octave/matlab if you want to use them.
Here is a code snipped that you could use to get the frequency of your columns in python
import numpy as np
omega = np.arange(2, 8).reshape(1, -1) # a row vector
t = np.linspace(0, 10, 1001).reshape(-1, 1) # a column vector
# this is a matrix as you described
matrix = np.angle(np.exp(1j * omega * t))
linearized = np.unwrap(matrix, axis=0);
# this function will get the frequency of your underlying signal
omega_ = np.mean(np.diff(linearized, axis=0), axis=0)

Then for the phase you could use
phase = np.mean(linearized - omega_ * t)