Удаление неразрешимых уравнений из недоопределенной системы

Моя программа пытается решить систему линейных уравнений. Для этого он собирает матрицу coeff_matrix и вектор value_vector и использует Eigen для их решения следующим образом:

Eigen::VectorXd sol_vector = coeff_matrix
        .colPivHouseholderQr().solve(value_vector);

Проблема в том, что система может быть как переопределена, так и недоопределена. В первом случае Эйген дает либо правильное, либо неправильное решение, и я проверяю решение с помощью coeff_matrix * sol_vector - value_vector.

Однако, пожалуйста, рассмотрите следующую систему уравнений:

a + b - c     =  0
        c - d =  0
        c     = 11
      - c + d =  0

В этом конкретном случае Эйген правильно решает три последних уравнения, но также дает решения для a и b.

Чего я хотел бы добиться, так это того, чтобы решались только те уравнения, которые имеют только одно решение, а остальные (здесь первое уравнение) оставались бы в системе.

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

Не могли бы вы предложить какой-нибудь хороший способ добиться этого?

Редактировать: обратите внимание, что в большинстве случаев матрица не будет квадратной. Я добавил сюда еще одну строку просто для того, чтобы отметить, что переопределение тоже может иметь место.


c++ eigen linear-algebra
person Michał Górny    schedule 10.07.2012    source источник

Ответы (3)


arrow_upward
6
arrow_downward

Я думаю, вам нужна разложение по сингулярным значениям (SVD), которая даст вам именно то, что вы хотите. После СВД «будут решаться уравнения, имеющие только одно решение», а решение псевдообратное. Это также даст вам нулевое пространство (откуда берутся бесконечные решения) и левое нулевое пространство (откуда возникает несоответствие, то есть нет решения).

person chaohuang    schedule 10.07.2012
comment
Не могли бы вы указать мне на какую-нибудь страницу или документ, который помог бы мне понять, как его использовать? Википедия много говорит о математике... - person Michał Górny; 10.07.2012
comment
На самом деле для изучения SVD требуется совсем немного математики, что является кульминацией линейной алгебры. Я только что погуглил несколько вводных, см. здесь и здесь. ХТН. - person chaohuang; 10.07.2012
comment
Спасибо. Я немного подвигался по пустым местам и собрал что-то, что я вставил в качестве решения (поскольку комментарии не годятся для реального кода). Не уверен, что я изобрел что-то реальное или просто случайный паттерн. - person Michał Górny; 10.07.2012

arrow_upward
2
arrow_downward

Основываясь на комментарии SVD, я смог сделать что-то вроде этого:

Eigen::FullPivLU<Eigen::MatrixXd> lu = coeff_matrix.fullPivLu();

Eigen::VectorXd sol_vector = lu.solve(value_vector);
Eigen::VectorXd null_vector = lu.kernel().rowwise().sum();

AFAICS, строки null_vector, соответствующие одиночным решениям, равны 0s, а строки, соответствующие неопределенным решениям, - 1s. Я могу воспроизвести это во всех моих примерах с порогом по умолчанию, который есть у Эйгена.

Однако я не уверен, что делаю что-то правильно или просто заметил случайную закономерность.

person Michał Górny    schedule 10.07.2012
comment
Если я правильно помню свою линейную алгебру, она кажется правильной. Я использовал подобные уловки на экзаменах по прикладной линейной алгебре, где нам разрешалось использовать Matlab, но приходилось отрабатывать решения шаг за шагом. - person Kuba hasn't forgotten Monica; 29.07.2012

arrow_upward
1
arrow_downward

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

Вот ссылка на то, что такое определитель и как его вычислить: http://en.wikipedia.org/wiki/Determinant

Обратите внимание, что исключение Гаусса также должно работать: http://en.wikipedia.org/wiki/Gaussian_elimination С помощью этого метода вы получите строки из 0, если существует бесконечное количество решений.

Изменить

Если матрица не квадратная, сначала нужно извлечь квадратную матрицу. Есть два случая:

  1. У вас больше переменных, чем уравнений: тогда у вас либо нет решения, либо их бесконечное количество.
  2. У вас больше уравнений, чем переменных: в этом случае найдите квадратную подматрицу ненулевого определителя. Решите эту матрицу и проверьте решение. Если решение не подходит, значит, у вас нет решения. Если решение подходит, это означает, что дополнительные уравнения линейно зависели от извлеченных.

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

Что касается исключения Гаусса, то оно должно работать напрямую с неквадратными матрицами. Однако на этот раз вы должны проверить, что количество непустых строк (т.е. строк с некоторыми ненулевыми значениями) равно количеству переменных. Если меньше, то у вас бесконечное количество решений, а если больше, то у вас нет решений.

person PierreBdR    schedule 10.07.2012
comment
Определитель действителен только в квадратных матрицах. - person Michał Górny; 10.07.2012
comment
@MichałGórny Я отредактировал свое решение для неквадратного случая. Но только квадратная матрица допускает единственное решение. Кроме того, ваш пример был с квадратной матрицей. - person PierreBdR; 10.07.2012
comment
Мой пример также представляет собой переопределенную матрицу, в которой две переменные имеют единственное решение. Все не так просто... Кстати, я думал о чем-то подобном, но боюсь, что все эти расчеты сделают решение ужасно неэффективным. - person Michał Górny; 10.07.2012