Реализация алгоритма дискриминанта для сокращенных многомерных признаков на Python

Вопрос или проблема

def algoritam2(G, I, F, max_iter, tol):
    """
    Алгоритм дискриминанта для уменьшенных многомерных признаков

    Параметры:
    G : numpy.ndarray
        Основная тензорная матрица признаков размерности J1 x J2 x ... x JN x K

    I : numpy.ndarray
        Целевые метки для основной матрицы G (K признаков в C классах)

    F : int
        Число дискриминационных признаков, которое нам нужно

    max_iter: int
        Максимальное число итераций при максимизации отношения следов

    tol: float
        Допуск для динамического критерия сходимости

    Возвращаемые значения:
    Psi : numpy.ndarray
        Матрица дискриминационной проекции, размерность L x F

    F_matrix : numpy.ndarray
        Матрица дискриминационных признаков, размерность F x K

    примечание: L = J1 * J2 * ... * JN
    """
    K = G.shape[-1] 
    N = G.ndim - 1  
    klase = np.unique(I)

    G_mean = np.mean(G, axis=-1)

    G_c_mean=np.zeros(G.shape[:-1]+(len(klase),))
    G_tilda=np.zeros(G.shape)

    for c in klase:
        ind=np.where(I==c)[0]
        G_c = G[..., ind]  # Образцы из класса c
        K_c = G_c.shape[-1]   # Число образцов в классе c
        G_c_mean[...,c] = np.mean(G_c, axis=-1)  # Среднее для класса c

        for i in ind:
            G_tilda[...,i]=G[...,i]-G_c_mean[...,c]  # Централизация класса c
        G_c_mean[...,c]=np.sqrt(K_c)*(G_c_mean[...,c]-G_mean)

    S_w = tl.base.unfold(G_tilda,N).T @ tl.base.unfold(G_tilda,N)
    S_b = tl.base.unfold(G_c_mean,N).T @ tl.base.unfold(G_c_mean,N)

    delta, psi = eigs(S_b,k=F,M=S_w,which="LM")

    br_iter=0
    uvjet=1

    while br_iter<max_iter and uvjet:
        psi_stari=psi
        br_iter+=1

        fi=np.trace(psi.T @ S_b @ psi)/np.trace(psi.T @ S_w @ psi) 
        delta, psi = eigs(S_b-fi*S_w,k=F,which="LM")
        delta, psi = eigs([email protected]@S_w@[email protected],k=F,which="LM")
        uvjet=(np.linalg.norm(psi-psi_stari)>=tol)

    F_matrix = psi.T @ tl.unfold(G, N).T
    return psi, F_matrix

Ответ или решение

Для того чтобы реализовать Алгоритм Дискриминанта для сокращенных многомерных признаков на языке Python, важно убедиться, что вы правильно понимаете теоретическую базу, лежащую в его основе, и используете соответствующие инструменты и библиотеки, такие как Tensorly, Scipy и Numpy.

Теория

Алгоритмы дискриминантного анализа часто применяются для выделения признаков в многомерных данных с целью выявления классов. Конкретно в вашем случае, алгоритм основывается на максимизации отношения следов двух ковариационных матриц: междуклассовой и внутриклассовой. Основная идея заключается в том, чтобы максимизировать межклассовое расстройство, одновременно минимизируя внутриклассовое, тем самым находя более четкие границы между различными классами в данных.

Пример

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

Применение

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

1. Использование Eigen-решения с комплексными значениями:
   Комплексные числа могут появляться вследствие применения функции eigs из Scipy, которая вычисляет только несколько собственных значений большой квадратной матрицы. Обратите внимание на использование параметра which="LM" – это может вызывать появление комплексных значений, если матрицы S_b и S_w становятся не совсем корректными.

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

3. Условие выхода из цикла:
   Проверьте логику выхода из цикла: (np.linalg.norm(psi-psi_stari) >= tol). Возможно, ваш порог tol слишком мал, что приводит к слишком жесткому критерию сходимости. Попробуйте увеличить этот порог или добавить дополнительные проверочные условия для более гибкой сходимости.

4. Проверка на ортогональность:
   Иногда полезно следить за ортогональностью получаемых векторов, и если она нарушается, применять процедуру ортогонализации, например с использованием метода Грама-Шмидта.

5. Нормализация данных:
   Убедитесь, что ваши данные должным образом нормализованы или стандартированы перед проведением анализа. Это может помочь избежать артефактов, связанных с численной масштабируемостью.

6. Обработка вывода:
   При многократных попарных вычислениях eigs может быть полезным дополнительно проверить результаты на наличие комплексной части или использовать операции, возвращающие только действительную часть в случае, если комплексная часть слишком мала.

Ниже представлен улучшенный подход с учётом описанных понятий:

import numpy as np
from scipy.sparse.linalg import eigs
import tensorly as tl

def algoritam2(G, I, F, max_iter, tol):
    K = G.shape[-1] 
    N = G.ndim - 1  
    klase = np.unique(I)

    G_mean = np.mean(G, axis=-1)

    G_c_mean = np.zeros(G.shape[:-1] + (len(klase),))
    G_tilda = np.zeros(G.shape)

    for c in klase:
        ind = np.where(I == c)[0]
        G_c = G[..., ind]
        K_c = G_c.shape[-1]
        G_c_mean[..., c] = np.mean(G_c, axis=-1)

        for i in ind:
            G_tilda[..., i] = G[..., i] - G_c_mean[..., c]

        G_c_mean[..., c] *= np.sqrt(K_c)
        G_c_mean[..., c] -= G_mean

    S_w = tl.unfold(G_tilda, N).T @ tl.unfold(G_tilda, N)
    S_b = tl.unfold(G_c_mean, N).T @ tl.unfold(G_c_mean, N)

    if np.linalg.matrix_rank(S_w) < S_w.shape[0]:
        S_w += np.eye(S_w.shape[0]) * tol  # Добавляем небольшую идентичную матрицу, если нужно

    delta, psi = eigs(S_b, k=F, M=S_w, which="LM")
    delta, psi = delta.real, psi.real

    for _ in range(max_iter):
        psi_stari = psi.copy()
        trace_quotient = np.trace(psi.T @ S_b @ psi) / np.trace(psi.T @ S_w @ psi)
        delta, psi = eigs(S_b - trace_quotient * S_w, k=F, which="LM")
        delta, psi = delta.real, psi.real
        if np.linalg.norm(psi - psi_stari) < tol:
            break

    F_matrix = psi.T @ tl.unfold(G, N).T
    return psi, F_matrix

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