Решение линейной системы уравнений с ограничениями на неизвестных
Я хотел бы решить систему линейных уравнений y=Uh для неизвестного вектора h, где у меня есть несколько ограничений на некоторые элементы $h$. Матрица U состоит из вектора u (длины N) и его задержанной версии.
Рассмотрим небольшой пример: y — вектор размера N, искаженный шумом, U — матрица размера $N \times 9$, где:
the 1st column of U is the vector u;
the 2nd column is one time step delayed vector u i.e. [0; u(1:end-1)];
the 3rd column is two time step delayed vector u i.e. [0; 0; u(1:end-2)];
the 4th column is a square of individual elements of 1st column;
the 5th column is a square of individual elements of 2nd column;
the 6th column is a square of individual elements of the 3rd column;
the 7th column is the product of individual elements of columns 1 and 2;
the 8th column is the product of individual elements of columns 1 and 3;
the 9th column is the product of individual elements of columns 2 and 3;
а $h$ — вектор размера $9$, где $N>9$. Здесь у меня мало ограничений на вектор $h$. Если первые три элемента — $a, b, c$, то следующие три элемента должны быть $pa^2, pb^2, pc^2$, где $p$ — ненулевой коэффициент масштабирования. Последними тремя элементами должны быть $2pab, 2pac, 2pbc$. Как использовать такие ограничения для параметров, используя другой параметр неизвестного вектора? Может ли кто-нибудь помочь мне с этим?
Решение h = pinv(U)*y не гарантирует желаемую структуру h. Как гарантировать, что решение имеет такую структуру и наилучшим образом аппроксимирует истинный шумный выходной сигнал y? Здесь неизвестными являются a,b,c,p.
p задано или неизвестно?
Система слишком детерминирована. Насколько вы уверены, что существует точное решение? (Вероятно, нет.) Вы пытаетесь выполнить аппроксимацию методом наименьших квадратов?
Ответы 1
Это не линейная система; это нелинейная система. Кроме того, вам почти наверняка придется принять более слабое определение решения, поскольку система переопределена и поэтому ее можно решить только методом наименьших квадратов. Неизвестных всего три (a, b, c), а уравнений N>9. На примере
import numpy as np
from scipy.optimize import least_squares
def make_h(abcp: np.ndarray) -> np.ndarray:
abc, (p,) = np.split(abcp, (-1,))
a, b, c = abc
# Here, I have few constraints on the vector $h$.
return np.concatenate((
# If the first three elements are $a, b, c$,
abc,
# then the next three elements should be $pa^2, pb^2, pc^2$.
p*abc**2,
# The last three elements should then be $2pab, 2pac, 2pbc$.
(
2*p*a*b,
2*p*a*c,
2*p*b*c,
),
))
def residuals(
abcp: np.ndarray, U: np.ndarray, y: np.ndarray,
) -> np.ndarray:
h = make_h(abcp)
# I would like to solve a system of linear equations $y=Uh$ for an unknown vector $h$.
yhat = U @ h
return yhat - y
def solve(
U: np.ndarray, y: np.ndarray,
) -> tuple[np.ndarray, np.ndarray]:
M = U.shape[1]//3
N = U.shape[0]
result = least_squares(
fun=residuals,
args=(U, y),
x0=np.ones(M + 1),
)
if not result.success:
raise ValueError(result.message)
return make_h(result.x), result.x
def main() -> None:
rand = np.random.default_rng(seed=0)
M = 3
N = 20 # $N>9$.
u = rand.random(size=N) # a vector u (length N)
padded = np.concatenate((np.zeros(M - 1), u))
delay = np.lib.stride_tricks.sliding_window_view(
x=padded,
window_shape=N,
)[::-1, :]
ua, ub, uc = delay
U = np.vstack((
delay, delay**2,
ua*ub, ua*uc, ub*uc,
)).T
# Make some fake data that stand a chance of being solvable
abcp_hidden = rand.random(size=M + 1)
h_hidden = make_h(abcp_hidden)
noise = rand.normal(size=N, scale=0.05)
# $y$ is a vector of size $N$
y = noise + U@h_hidden
h, (*abc, p) = solve(U, y)
print('h =', h)
print('abc =', abc)
print('p =', p)
if __name__ == '__main__':
main()
h = [2.70105108e-02 1.47198194e-01 6.76129121e-01 4.22152723e-04
1.25374428e-02 2.64522905e-01 4.60118054e-03 2.11347169e-02
1.15177095e-01]
abc = [0.02701051078001584, 0.14719819396169204, 0.6761291207038694]
p = 0.5786340691198476
Благодарю за ваш ответ. Я отредактировал свой вопрос для большей ясности. Матрица U не случайна, а имеет определенную структуру. Фактически, он создан из вектора u и его задержанных версий. Масштабный коэффициент p здесь также неизвестен. Да, система переопределена, однако я ищу решение, которое лучше всего приближает y и при этом имеет желаемую структуру h.
Я знаю, что U не случайно, но вы (до сих пор) не смогли предоставить минимально воспроизводимый пример, поэтому мне пришлось кое-что заполнить. Минимальный воспроизводимый пример должен иметь возможность запускаться от начала до конца после копирования и вставки.
Тем не менее: от скуки я отредактирую этот ответ, чтобы продемонстрировать, как будет выглядеть U-образная конструкция.
У меня есть измерения y и u в системе, где выходной сигнал y зашумлен. Эти измерения поступают из системы, которая в идеале должна следовать описанной мной математической модели, т.е. на временном шаге k выход y_k и вход u_k связаны следующим образом: y_k = au_k + bu_{k-1} + cu_{k-2 } + p(au_k + bu_{k-1} + cu_{k-2} )^2. Я пишу это уравнение как линейную систему уравнений y = Uh для k = 1,...,N. и мне нужно найти h. На данный момент решение h = pinv(U)*y не обеспечивает требуемую структуру h. Как мне обеспечить это?
Здесь вы можете предположить, что u — случайный вектор, который используется для построения матрицы U, а y следует математической модели с заданными a, b, c и p. Теперь к y добавляется шум, и я хотел бы найти заданные a, b, c и p, используя шумовое измерение y и матрицу U.
Спасибо за исправление в вашем решении. Вектор решения h отклоняется от желаемой структуры с увеличением шума. Это связано с тем, что функция «решить» использует метод наименьших квадратов и не ограничивает требуемую структуру решения. Кроме того, предположим, что я фиксирую N = 1000 и выбираю первые 500 выборок для создания y как: y(1:500) = U(1:500,:)*h + Noise1. Затем я решаю для h. Используя этот h, я вычисляю y(501:1000), используя U(501:1000), тогда возникает значительное отклонение. Обычно это набор обучающих и тестовых данных.
Это очевидно. Единственный раз, когда можно надеяться на точное решение, — это точно определенный или недоопределенный случай. Для реалистичных данных ваше решение гарантированно будет неидеальным, кроме как в смысле наименьших квадратов.