Есть ли способ генерировать случайные решения неквадратных линейных уравнений, желательно на питоне?

Система линейных уравнений с большим количеством переменных, чем уравнений, известна как недоопределенная система.

Недоопределенная линейная система либо не имеет решений, либо имеет бесконечно много решений.

...

Существуют алгоритмы, позволяющие определить, есть ли решения у недоопределенной системы, и, если они есть, выразить все решения в виде линейных функций от k переменных (таких же k, как указано выше). Самый простой из них — исключение Гаусса.

Как вы говорите, многие функции, доступные в библиотеках (например, np.linalg.solve), требуют квадратной матрицы (то есть n уравнений для n неизвестных), то, что вы ищете, — это реализация Исключение Гаусса для неквадратных линейных систем.

Это не «случайно», а np.linalg.lstsq (наименьший квадрат) будет решать неквадратные матрицы:

Верните решение по методу наименьших квадратов линейному матричному уравнению.

Решает уравнение a x = b путем вычисления вектора x, минимизирующего евклидову 2-норму || б - а х ||^2. Уравнение может быть недоопределено, хорошо или переопределено (т. е. количество линейно независимых строк может быть меньше, равно или больше числа его линейно независимых столбцов). Если a квадратное и имеет полный ранг, то x (за исключением ошибки округления) является «точным» решением уравнения.

Для получения дополнительной информации см.: решение Ax = b для неквадратной матрицы A с использованием python

Поскольку у вас есть недоопределенная система уравнений (слишком мало ограничений для ваших решений или меньше уравнений, чем переменных), вы можете просто выбрать некоторые произвольные значения для x3 и x4 и решить систему в x1, x2 (это имеет 2 переменные/2 уравнения).

Вам просто нужно будет проверить, что получившаяся система не является противоречивой (т. е. не допускает решения) и что нет повторяющихся решений.

Например, вы можете зафиксировать x3=0 и выбрать случайные значения x4 для создания решений ваших уравнений в x1, x2.

Вот пример генерации 10 "случайных" решений

n = 10
x3 = 0
X = []
for x4 in np.random.choice(1000, n):
  b = np.array([[6-7*x3+2*x4],[-2+2*x3-5*x4]])
  x = np.linalg.solve(a, b)
  X.append(np.append(x,[x3,x4]))

# check solution nr. 3
[x1, x2, x3, x4] = X[3]
3*x1 - 2* x2 + 7*x3 - 2*x4
# output:  6.0
1*x1 + 3* x2 - 2*x3 + 5*x4
# output: -2.0

Спасибо за ответы, которые помогли мне и указали мне правильное направление.

Теперь у меня есть простое пошаговое решение моей проблемы для произвольного k<n.

1. Найдите одно решение всех приведенных уравнений. Это можно сделать с помощью

 solution_vec = numpy.linalg.lstsq(A,b)

это дает решение, как показано в ответе ukemis. В моем примере выше Матрица A равна коэффициентам уравнений в левой части, b представляет вектор в правой части.

2. Определите нулевое пространство вашей матрицы A.

Это все векторы v такие, что скалярное произведение v*A_i = 0 для каждой(!) строки A_i из A. Следующая функция найдена в этом потоке можно использовать для получения представителей нулевого пространства A:

def nullSpaceOfMatrix(A, eps=1e-15):
    u, s, vh = scipy.linalg.svd(A)
    null_mask = (s <= eps)
    null_space = scipy.compress(null_mask, vh, axis=0)
    return scipy.transpose(null_space)  

3. Генерируйте столько (N) "случайных" линейных комбинаций (то есть со случайными коэффициентами) solution_vec и результирующих векторов нулевого пространства матрицы, сколько хотите! Это работает, потому что скалярное произведение является аддитивным, а векторы нулевого пространства имеют скалярное произведение 0 на векторы уравнений. Эти линейные комбинации всегда должны содержать solution_vec, например:

linear_combination = solution_vec + a*null_spacevec_1 + b*nullspacevec_2...

где a и b могут быть выбраны случайным образом.