Принадлежит ли точка конечно порождённому конусу?

[B@R5uk]

У меня имеется N ненулевых векторов a с индексом в пространстве размерности тоже N:

$$a_k\in\mathbb R^N,\qquad\qquad 1\le k\le N$$

Пусть для начала они все линейно независимые (чтобы не возиться с вырожденными случаями, пока задача не решена), то есть образуют базис. Если рассмотреть множество точек, являющихся концами вектора x линейной комбинации этих исходных векторов с неотрицательными коэффициентами, то получится конечно порождённый конус K с вершиной в начале координат:

$$K=\left\langle\left{}\;x=\sum\limits_{k=1}^{N}\lambda_ka_k\;\right|\;\lambda_k\ge 0,\;1\le k\le N\;\right\rangle$$

Теперь, у меня также имеется вектор b:

$$b\in\mathbb R^N$$

Требуется выяснить: принадлежит ли точка, задаваемая этим вектором b конусу K?

Я знаю, что существует лемма Фаркаша об альтернативе, которая как раз рассматривает этот вопрос в самом общем случае. Однако, она предоставляет критерии для проверки, типа "если да, то существует то-то, если нет, то существует кое-что другое", но не говорит как найти эти "то" и "кое-что", являющиеся свидетелями того или другого случая. Для поиска придётся воспользоваться линейным программированием. Но решать задачу ЛП для моего примитивного случая, когда количество базовых векторов равно (или даже меньше, если они ЛЗ) размерности пространства, кажется как-то уж через чур. Я подозреваю, что в этом простом виде задачу можно решить методами базовой линейной алгебры и/или аналитической геометрии.

Например, в случае плоскости, всё совсем элементарно. Надо проверить, что угол, который составляет вектор b с базовыми векторами a меньше (или равен) углу между этими векторами:

$$\begin{cases}\lbig(a_1,\;b\rbig)\ge\lbig(a_1,\;a_2\rbig)\\\lbig(a_2,\;b\rbig)\ge\lbig(a_1,\;a_2\rbig)\end{cases}$$

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

Случай трёхмерного пространства заметно сложнее. Первое, что приходит на ум — это найти вектор c, который лежит строго внутри конуса, например, взяв сумму базовых векторов a с весом 1:

$$c=a_1+a_2+a_3$$

А затем проверить, что вектора b и c лежат по одну сторону от каждой плоскости трёхгранного угла, ограничивающего конус K:

$$\begin{cases}\lbig(b,\;n_1\rbig)\lbig(c,\;n_1\rbig)\ge 0\\\lbig(b,\;n_2\rbig)\lbig(c,\;n_2\rbig)\ge 0\\\lbig(b,\;n_3\rbig)\lbig(c,\;n_3\rbig)\ge 0\end{cases}$$

Здесь n с индексом — это вектора-нормали к каждой из плоскостей, которые можно вычислить, например, взяв векторное произведение двух базовых векторов a с соответствующими индексами. Вектора в этих неравенствах можно не нормировать, так для проверки важен только знак.

Этот подход кажется обобщаемым на случай произвольной размерности, однако не совсем понятно, как именно это делать, особенно в плане вычисления векторного произведения. Вернее, тут надо искать нормаль к каждой гиперграни гипергранного угла, и я не понимаю, как это делать.

Задача является вычислительно-прикладной с высокой степенью задействованности (читай: крутится внутри самого внутреннего цикла), поэтому очень важно, чтобы она считалась быстро. Подскажите, пожалуйста, есть ли какие-то более оптимальные в вычислительном плане подходы (ну, или хотя бы просто другие подходы)?

[vpb]

Какое есть растение семейства портулаковых, 8 букв, первая "п", последняя "к" ? Портулак. Так и тут: разложить по базису, и посмотреть, положительны ли коэффициенты. А проще, думаю, нельзя. Для разложения же по базису надо решить систему линейных уравнений. Для чего есть методы (гауссово исключение, причем с выбором ведущего элемента, для численной устойчивости). Если же векторов больше, чем размерность, то всё сложно. Надо линейное программирование включать, но в нём я не разбираюсь.

[B@R5uk]

vpb в сообщении #1685195 писал(а):

...разложить по базису и посмотреть, положительны ли коэффициенты.

Точно! Что-то я затупил сильно. Конус является главным ортантом базиса векторов, очевидно, что любая точка внутри главного ортанта этого базиса будет иметь все положительные координаты. И вырожденные случаи легко решаются сравнением ранга обычной и расширенной матриц. Спасибо за помощь!

[B@R5uk]

B@R5uk в сообщении #1685197 писал(а):

И вырожденные случаи легко решаются сравнением ранга обычной и расширенной матриц.

На самом деле тут есть нюанс. Если вектор b лежит вне линейной оболочки (ЛО) базисных векторов a, то он, понятное дело, будет лежать и вне конечно порождённого конуса K. Если же он лежит в ЛО, то тут ситуация становится сложнее: мы имеет задачу общего вида, когда количество векторов, задающих конус, больше размерности пространства (которым является подпространство ЛО). Пока не знаю, нужно ли мне будет рассматривать такой случай.

В случае же, когда размерность пространства больше числа векторов, нужно искать псевдообратную матрицу, что можно легко сделать, например, через сингулярное разложение. Через него же считается и ранг матрицы (то есть проверяется линейна независимость векторов). Тестовый код поиграться в Matlab'е:

Используется синтаксис Matlab M

clc
clearvars
format compact
%format long
format short

%   Размерность пространства
N = 5;
%   Количество векторов базиса
M = 3;
%   Количество изучаемых векторов
K = 13;

%   Матрица стобцов координат базиса
a = randn (N, M);
%   Экономное сингулярное разложение базисной матрицы
%   u s v' = a
[u s v] = svd (a, 'econ');
%   Если базис вырожден, то либо изучаемый вектор будет лежать вне
%   линейной оболочки базиса и, поэтому, вне образуемого конечно
%   порождённого конуса; либо он будет лежать в этой оболочке и задача
%   приобретёт общий характер задачи ЛП.
if s (end) / s (1) < eps * N
    error ('Basis is degenerate.')
end
%   Изучаемые вектора все лежат в подпространстве базистных векторов
b = a * randn (M, K);
%   Координаты изучаемых векторов в базисе
%c = a \ b;
c = v * s ^ -1 * u' * b;
%   Проверка, что разложение полное (невязка,
%   заметно большая машинной точности, будет свидетельсвовать об обратном)
disp (ceil (abs (sum ((a * c) .^ 2) - sum (b .^ 2)) / eps))
%   Строка с единицами для векторов в конусе
flags = prod (double (0 <= c));
disp (flags)
%   Исходные координаты векторов в конусе
disp (b (:, 0 ~= flags))
disp (' ')
%   Координаты разложения векторов в конусе по базису
disp (c (:, 0 ~= flags))