Введение
Мы знаем, что для системы линейных уравнений, в которой количество уравнений меньше числа переменных, должно быть бесконечное число решений. И что количество свободных переменных должно быть равно количеству переменных минус количество уравнений.
Но как мы можем использовать этот факт? Например, для задачи оптимизации с ограничениями линейного равенства можно было бы надеяться использовать ограничения для уменьшения количества переменных. А затем иметь возможность использовать этот новый сокращенный набор переменных в алгоритме неограниченной оптимизации.
Как мы можем этого добиться? Введите нулевое пространство.
Нулевое пространство
Давайте определим нашу недоопределенную систему из M уравнений и N переменных как:
Нулевое пространство матрицы A определяется как Z. Это набор N-M векторов-столбцов, которые «обнуляют» A. Каждый из этих независимых векторов-столбцов дает нулевой вектор при умножении на A слева.
Поскольку N-M векторов-столбцов в Z дают нулевой вектор при умножении на A, линейная комбинация этих векторов также будет делать то же самое.
Таким образом, мы сократили количество переменных с N до N-M переменных в векторе u.
Однако это справедливо только для однородной линейной системы. Это системы с нулевой правой частью. Но, имея единственное решение нашей линейной системы, мы можем преобразовать ее в однородную систему:
И мы можем сделать вывод:
Нулевое пространство в SciPy
В SciPy есть функция, которая вычисляет ортонормированный базис нулевого пространства Z с использованием разложения по сингулярным значениям. Он автоматически определит количество столбцов в Z на основе вычисленных сингулярных значений. Z будет иметь столбцы N-ранга (A). Если все строки линейно независимы, то Rank(A)=M.
import numpy as np from scipy.linalg import null_space np.random.seed(0) # Our linear system M,N = 100,200 A,b = np.random.randn(M,N),np.random.randn(M) # Calculate Z Z = null_space(A) # An easy way to calculate a single solution to Ax = b: # This method will give the solution with the minimal norm ||x||. x0 = np.linalg.pinv(A)@b # Function to calculate x(N) from u(N-rank(A)) x_from_u = lambda u: Z@u+x0 # Generate a new x from a random u xnew = x_from_u(np.random.randn(Z.shape[1])) # Check that the linear system is still solved. assert(np.allclose(A@xnew,b))
Обсуждение
Хотя SVD является наиболее точным способом численного расчета базиса нулевого пространства, может возникнуть некоторая ошибка. Мы можем убедиться в этом, сравнив A.dot(Z) с 0. В приведенном выше коде он имеет порядок 1e-14. Но оно будет расти вместе с величиной тебя. Если u станет слишком большим, то Ax=b больше не будет в достаточной степени удовлетворяться.
Таким образом, рекомендуется убедиться, что решение x0 близко к области x, которую вы исследуете. Следующий код находит решение x0 вблизи указанной точки v.
from scipy.linalg import solve
# Minimize ||x-v||**2 subject to Ax = b. Using lagrange multipliers.
# I.e. finds a solution to Ax=b close to v.
def get_solution(v,A,b):
M,N=A.shape[0],A.shape[1]
mat = np.block([[np.eye(N),A.T],[A,np.zeros((M,M))]])
y = np.concatenate([v,b])
ret = solve(mat,y,assume_a="sym")[0:N]
return ret
x0 = get_solution(np.ones(N),A,b)
В заключение
Функцию null_space() в SciPy можно использовать для перехода к новому, меньшему набору переменных. Если между вашими переменными существует связь, описываемая линейной системой.
Помимо прочего, его можно использовать для реализации ограничений линейного равенства без привязки к алгоритму оптимизации с ограничениями.
Это не тот метод, который часто используется, и следует опасаться числовых неточностей. Т.е. часто проверяйте, что линейная система удовлетворяется (||A(x0+Zu)-b||**2 мала).
Дополнительный пример с ограниченной линейной регрессией
В качестве примера давайте решим следующую задачу оптимизации двумя способами. Сначала воспользуемся стандартным методом множителей Лагранжа. Затем, изменив переменные с помощью null_space().
Используя множители Лагранжа, каждое ограничение вводит новую переменную, которую необходимо решить. Полученная линейная система также больше не является положительно-определенной.
Изменяя переменные с использованием пустого пространства, мы удаляем по одной переменной для каждого ограничения. И полученная линейная система будет положительно определенной. Но, как упоминалось выше, могут возникнуть проблемы с небольшими числовыми неточностями в вычислении null_space(), которые растут по мере удаления от нашего единственного решения x0.
from scipy.linalg import null_space,solve import numpy as np ## ## Generate some data #### N,M,NC=200,50,20 # Number of datapoints, variables and constraints. X=np.concatenate([np.random.randn(N,M-1),np.ones((N,1))],axis=1) [email protected](M) +np.random.randn(N) # Generate some constraints A,b = np.random.randn(NC,M),np.random.randn(NC) # Let's also add some l2-regularization... l2_reg = np.eye(M) l2_reg[-1,-1] = 0.0 # Only on weights, not on the bias. # Define Q and c in the loss function. # This is from the standard form of quadratic programming loss function. Q=X.T@X+l2_reg c = -X.T@y ## ## Solution using Lagrange multipliers #### mat = np.block([[Q,A.T],[A,np.zeros((NC,NC))]]) vec = np.concatenate([-c, b]) sol1 = solve(mat,vec,assume_a="sym")[0:M] ## ## Solution using null space to change variables. #### Z = null_space(A) x0 = np.linalg.pinv(A)@b #Just plug in our formula x=Zu+x0 in the loss function. And get: Q_new = Z.T@Q@Z c_new = Z.T@(c+Q@x0) u = solve(Q_new,-c_new,assume_a="pos") sol2 = Z@u+x0 ## ## Check the two solutions #### assert(np.allclose(sol1,sol2))
