Геометрический паспорт матрицы: извлекаем инварианты из n-мерных симплексов / Песочница / Хабр

В чем идея?

Традиционный спектральный анализ характеризует оператор в целом, но не всегда наглядно отражает геометрию. Предложенный подход позволяет перейти к дескрипции на основе трех фундаментальных векторов: аналитического отклика $\vec{s}$ , метрического радиус-вектора $\vec{r}$ и вектора асимметрии масс $\vec{p}$ .

Мы должны получить дескриптор , содержащий глобальные инварианты :
- радиус описанной вокруг геометрического тела гиперсферы.
$G \in [0,1[$ - мера "кучности" расположения точек симплекса на гиперсфере или глобальная асимметрия симплекса.
и набор локальных инвариантов $(\alpha_i, \beta_i, g_i)$ для каждой -вершины симплекса, где:
— количество вершин симплекса равно .
$\alpha_i$ — угол "раскрытия вершины".
$\beta_i$ — радиальный эксцентриситет (угол смещения центра).
— вклад вершины в общую асимметрию.

Алгоритм вычисления для i-компонента дескриптора

Алгоритм состоит из 13 шагов. Он базируется на инверсии исходной матрицы и анализе метрических характеристик получившегося симплекса (тела).
Входные данные:
- Матрица , полученная сменой базиса для -вершины, при исходная матрица ,
в алгоритме опущен нижний индекс.
- Вектор $\vec{e}$ , все координаты которого равны 1 : $e_k = 1, \text{ где } k=1,2..n$ .
- Вектор $\vec{q}$ , координаты которого равны сумме квадратов элементов строк матрицы : $q_k = \sum^n_{j=1} a_{kj}^2 ,\text{ где } k=1,2..n.$

Шаги алгоритма :

$\vec{s} = A^{-1}\vec{e}$
$S = |\vec{s}|$
$\vec{d} = A^{-1}\vec{q}$
$\vec{r}=-\vec{d}/2$
$R = |\vec{r}|$
$\cos(\beta) = (\vec{s} \cdot \vec{r}) / (S \cdot R)$
$\vec{s}\;' = \frac{1}{S^2}\vec{s}$
$L = |\vec{s}\;'|$
$\cos(\alpha) = \frac{L}{R} + \cos(\beta)$
$\vec{t} = A^T \vec{e} + \vec{r}$
$\vec{p} = (-(n+1)\vec{r} - \vec{t}) / (n+1)$
$G = |\vec{p}|/R$
$g = \vec{p} \;\vec{r}/R^2$

Выходные данные:

- — радиус гиперсферы, глобальный инвариант дескриптора
- $G \in [0,1[$ - мера "кучности" точек симплекса на поверхности гиперсферы, глобальный инвариант дескриптора. После получения всех n+1 компонентов дескриптора должно выполняться равенство

$\sum^n_{i=0} g_{i} = (n+1)G^2$ .

- $\alpha \in [0,\pi[, \beta \in [0,\pi] , \alpha \ne \beta$ — условия невырожденности матрицы для каждого -компонента дескриптора.
$R,G,\alpha,\beta,g$ , вычисленные для каждого -компонента есть инварианты относительно операций в текущей матрице :
- перестановка строк или столбцов,
- умножение на матрицу поворота.
Угол $\alpha$ дополнительно инвариантен относительно операции умножения строки на положительное число.

Геометрическая интерпретация

На Рис. 1 продемонстрирована связь векторов алгоритма для A(2,2) , имеющей векторы-строки $\vec{a_1}, \vec{a_2}$ . Вектор $\vec{r}$ определяет положение центра описанной сферы относительно текущего базиса , Вектор $\vec{s}$ - аналитический отклик системы на единичное воздействие .Угол $\alpha$ - для размерности 2 равен углу между векторами $\vec{a_1}, \vec{a_2}$ и отражает структурную целостность («раскрытие» вершин) . Из Рис.1 очевидно , что угол $\alpha$ не изменит своего значения при умножении $\vec{a_1}$ и $\vec{a_2}$ на положительное число, а угол $\beta$ — радиальный эксцентриситет.

Тест № 1 Несимметричный тетраэдр

Рассмотрим несимметричный тетраэдр, заданный матрицей .

= $\begin{pmatrix} 1 & 0.2 & 0.1 \\ 0.2 & 1 & -0.1 \\ -0.1 & 0.1 & 1.5 \end{pmatrix}$

Тест подтверждает сохранение глобальных параметров при переносе начала координат в любую из 4-х вершин.

Верификация . Сумма $\sum^n_{i=0} g_i = 0.7691$ в точности равна теоретическому значению .

Ниже приведен код для Теста №1 на Python . Этот код не оптимизирован и служит только проверки логической корректности алгоритма.

import numpy as np

def get_matrix_passport(A):
    n_dim = A.shape[0]
    e = np.ones(n_dim)
    q = np.sum(A**2, axis=1)
    
    A_inv = np.linalg.inv(A)
    
    s = A_inv @ e
    S_norm = np.linalg.norm(s)
    
    d = A_inv @ q
    r = -d / 2
    R = np.linalg.norm(r)
    
    cos_beta = np.dot(s, r) / (S_norm * R)
    beta_deg = np.degrees(np.arccos(np.clip(cos_beta, -1, 1)))
    
    s_prime = s / (S_norm**2)
    L = np.linalg.norm(s_prime)
    cos_alpha = (L / R) + cos_beta
    # Safety clip for arccos
    alpha_val = np.clip(cos_alpha, -1, 1)
    alpha_deg = np.degrees(np.arccos(alpha_val))
    
    # t = сумма строк + r
    t = np.sum(A, axis=0) + r
    p = (-(n_dim + 1) * r - t) / (n_dim + 1)
    G = np.linalg.norm(p) / R
    
    g = np.dot(p, r) / (R**2)
    
    return R, G, alpha_deg, beta_deg, g

def shift_basis(A, i):
    if i == 0:
        return A
    target_vertex = A[i-1].copy()
    new_A = A - target_vertex
    all_vertices = np.vstack([-target_vertex, new_A])
    final_vertices = np.delete(all_vertices, i, axis=0)
    return final_vertices

A_initial = np.array([
    [1.0, 0.2, 0.1],
    [0.2, 1.0, -0.1],
    [-0.1, 0.1, 1.5]
])

results = []
for i in range(4):
    A_current = shift_basis(A_initial, i)
    res = get_matrix_passport(A_current)
    results.append((f"V{i}", *res))

print(f"{'Вершина':<8} | {'R':<8} | {'G':<8} | {'alpha':<8} | {'beta':<8} | {'g_i':<8}")
print("-" * 60)
for r in results:
    print(f"{r[0]:<8} | {r[1]:.4f} | {r[2]:.4f} | {r[3]:6.2f}° | {r[4]:6.2f}° | {r[5]:.4f}")

Тест № 2 Матрица Гильберта H(8,8)

Метод протестирован на экстремально плохо обусловленной матрице Гильберта размерностью 8.
Несмотря на близость к машинному нулю, алгоритм сохраняет стабильность инвариантов
и , что подтверждает его высокую численную устойчивость по сравнению со стандартными спектральными методами.

Верификация. Фактическая сумма $\sum g_i = 5.05600$ , теоретическое значение .

Ниже приведен код для Теста №2 на Python . Этот код не оптимизирован и служит только проверки логической корректности алгоритма.

import numpy as np

# Используем функцию расчета дескриптора из предыдущего примера
# (get_matrix_passport и shift_basis должны быть определены выше)

def test_hilbert_8x8():
    n = 8
    # Генерируем матрицу Гильберта 8x8
    H = np.fromfunction(lambda i, j: 1.0 / (i + j + 1), (n, n))
    
    print(f"Тест №2: Матрица Гильберта {n}x{n}")
    print(f"{'Вершина':<10} | {'R':<10} | {'G':<10} | {'alpha':<10} | {'beta':<10}")
    print("-" * 65)
    
    results = []
    for i in range(n + 1):
        A_curr = shift_basis(H, i)
        res = get_matrix_passport(A_curr)
        if res:
            results.append(res)
            # Выведем только некоторые вершины для экономии места
            if i in [0, 1, 7, 8]:
                r, g_glob, a, b, gi = res
                print(f"V{i:<9} | {r:.6f} | {g_glob:.6f} | {a:8.4f}° | {b:8.4f}°")

    # Проверка контрольной суммы
    sum_gi = sum(res[4] for res in results)
    theory_val = (n + 1) * (results[0][1]**2)
    print(f"\nКонтроль: Сумма g_i = {sum_gi:.5f}, (n+1)G^2 = {theory_val:.5f}")

test_hilbert_8x8()

Критерий идентичности

Два матричных дескриптора и признаются идентичными, если выполняются следующие условия:

Совпадение глобальных инвариантов: и с заданной точностью $\epsilon$ .
Совпадение локальных характеристик: неупорядоченные множества троек локальных компонентов $\{\alpha, \beta, g\}_A$ и $\{\alpha, \beta, g\}_B$ эквивалентны.

Заключение

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

Инвариантен к вращению, перестановке строк и выбору базовой вершины.
Устойчив к ошибкам округления на плохо обусловленных матрицах.
Целостен: локальные вклады вершин математически связаны с глобальной асимметрией фигуры .

Предложенный подход может быть использован в задачах идентификации систем, классификации данных и анализе устойчивости многомерных структур. Параметр $G \in [0,1[$ , обладающий высокой устойчивостью при вычислениях, может быть использован в качестве компонента первичного ключа при поиске .