Разрешение многомерной кривой: введение

Это сообщение воспроизведено из моего блога в Nirpy Research.

Предположим, вы измеряете спектр выборки, содержащей несколько компонентов, и хотите выделить вклад в спектр каждого компонента. Разрешение многомерной кривой (MCR) — это группа методов, разработанных для решения этой проблемы.

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

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

В процессе мы применим концепции, которые уже были рассмотрены в этом блоге, такие как градиентный спуск и алгоритм Кеннард-Стоун, и изучим несколько новых.

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

Спектры синтетических смесей

Данные, использованные в этом посте, доступны в нашем репозитории Github. Спектральные смеси были сгенерированы из данных, доступных на сайте RRUFF Project, содержащих рамановские спектры большого количества минералов.

Мы выбрали три минерала: кварц, ортоклаз и альбит и написали скрипт, который генерирует набор спектров смесей со случайными концентрациями этих трех чистых компонентов. Единственным ограничением является то, что (нормализованная) сумма концентраций должна равняться единице. Скрипт для генерации синтетических спектров пока доступен только нашим замечательным сторонникам Patreon.

Сценарий генерирует файл csv (который доступен на Github) из 300 случайных смесей и добавляет в конце три чистых спектра, всего 303 спектра. Давайте посмотрим на это

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
 
# Load data
data = pd.read_csv("Raman_spectra_mixtures.csv")
 
# Read the concentrations arrays
conc = data[["Albite", "Quartz", "Orthoclase"]].values[:-3,:]
 
# Read spectra arrays of the mixtures
spectra = data.values[:-3,4:1004]
 
# read pure spectra
pure_spectra = data.values[-3:,4:1004]
 
# Read wavenumbers
wn = data.keys()[4:1004].astype('float32')
 
with plt.style.context(('seaborn-whitegrid')):
    fig, ax = plt.subplots(figsize=(8, 5))
    for i,j in enumerate(["Albite", "Quartz", "Orthoclase"]):
        plt.plot(wn, pure_spectra[i,:], lw=2, label=j)
plt.xlabel("Wavenumber (cm$^{-1}$)", fontsize=14)
plt.ylabel("Raman spectra", fontsize=14)
plt.legend(fontsize=14)
plt.show()

Три компонента имеют существенно разные спектры, но с некоторыми перекрывающимися характеристиками.

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

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

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

Нахождение чистых спектров по концентрациям

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

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

Градиентный спуск — это процесс уменьшения ошибки путем изменения угаданных чистых спектров в сторону уменьшения градиента. Количество изменений (то есть размер шага) определяется скоростью обучения. Давайте определим нашу модифицированную версию градиентного спуска в специальной функции.

def gradient_descent(learning_rate, spectra, concentrations, n_iter=1000, pure_spectra_guess=None):
 
    if pure_spectra_guess is None:
        pure_spectra_guess = np.random.random((concentrations.shape[1], spectra.shape[1]))
    elif pure_spectra_guess == 'zero':
        pure_spectra_guess = np.zeros((concentrations.shape[1], spectra.shape[1]))
    else:
        pass
 
    error = [] # error
    for l in range(n_iter):
        # Calculate estimated spectra
        est_spectra = np.dot(concentrations, pure_spectra_guess)
        # Calculate errors
        e = spectra - est_spectra
        error.append(e.sum())
        # print(l, e.sum())  # optional print
        # Calculate gradient vector
        g = -np.dot(conc.T,e)/spectra.shape[0]
        # Gradient descent -&gt; update the guess
        pure_spectra_guess -= learning_rate*g
        pure_spectra_guess = np.abs(pure_spectra_guess) #positivity constraint
 
    return(pure_spectra_guess, error)

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

Хорошо, давайте запустим эту функцию и построим график результата

eta = 0.1
pure_spectra_gd, error = \
    gradient_descent(eta, spectra, conc, n_iter=5000, pure_spectra_guess=None)
 
with plt.style.context(('seaborn-whitegrid')):
    fig, ax = plt.subplots(figsize=(8, 5))
    plt.plot(pure_spectra_gd[0,:], lw=4, label = 'Original pure spectrum')
    plt.plot(pure_spectra[0,:], 'r', label = 'Guessed spectrum')
 
    plt.xlabel("Wavenumber (cm$^{-1}$)", fontsize=14)
    plt.ylabel("Raman spectra", fontsize=14)
    plt.legend(fontsize=14)    
    
    plt.show()

При достаточном количестве итераций предполагаемые чистые спектры сходятся к правильной форме, что мы и проверили, сравнив их с реальными чистыми спектрами.

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

Нахождение концентраций по чистым спектрам

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

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

from sklearn.linear_model import LinearRegression
est_conc = np.zeros_like(conc)
for i in range(conc.shape[0]):
    lr = LinearRegression()
    lr.fit(pure_spectra.T, spectra[i,:])
    est_conc[i,:] = lr.coef_

Концентрации — это просто коэффициенты линейной регрессии, по одному для каждого образца. Вы можете убедиться, сравнив массивы, что расчетные концентрации идентичны или очень близки к фактической концентрации.

Разрешение многомерной кривой

Это оно. Если вы сделали это так далеко, вот награда. Оценим и концентрации, и чистые спектры, исходя из спектральных смесей.

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

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

К счастью для нас, есть удобная библиотека Python, которую мы можем использовать для решения этой проблемы, которая называется pyMCR. В оставшейся части этого раздела я опишу метод решения проблемы MCR с нашим игрушечным набором данных, используя одну из функций pyMCR. Нам нужно несколько импортов, которые мы опишем через минуту.

from pymcr.mcr import McrAR
from pymcr.regressors import NNLS
from pymcr.constraints import ConstraintNonneg, ConstraintNorm
import kennard_stone as ks

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

Проблема нахождения трех самых разных наборов данных — это проблема, которую мы уже решили, когда обсуждали алгоритм Кеннарда-Стоуна. Алгоритм KS специально разработан для поиска ряда выборок в наборе данных, которые находятся дальше всего друг от друга, для некоторого определения наиболее удаленного, то есть для некоторого определения метрики расстояния. Эти детали обсуждаются в посте, указанном выше, поэтому здесь мы просто применяем разделение Кеннарда-Стоуна к нашему набору данных, чтобы извлечь три спектра, которые нам нужны.

# Estimate pure spectra picking the three most different spectra in the set
est_pure_spectra, X_test = ks.train_test_split(spectra, train_size = 3)

Теперь воспользуемся функцией McrAR, следуя примерам, приведенным на этой странице.

mcrar = McrAR(max_iter=1000, st_regr='NNLS', c_regr='OLS',
              c_constraints=[ConstraintNonneg(), ConstraintNorm()])
mcrar.fit(spectra, ST=est_pure_spectra)

Функция требует оценки либо концентраций, либо чистых спектров. Мы используем чистые спектры, рассчитанные по алгоритму КС.

Кроме того, функция требует двух процедур регрессии для чистого спектра (st_regr) и концентрации (c_regr). Мы используем неотрицательный наименьший квадрат (NNLS) и обычный метод наименьших квадратов (OLS) соответственно (кстати, подход OLS идентичен тому, который мы объяснили выше для калибровки концентраций).

Наконец, мы указываем ограничения для концентраций, которыми в нашем случае являются неотрицательность (концентрации являются положительными величинами) и нормировка (сумма концентраций равна единице).

После этого подгоняем спектральные смеси.

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

with plt.style.context(('seaborn-whitegrid')):
    fig, ax = plt.subplots(figsize=(8, 5))
    plt.plot(wn, mcrar.ST_.T, 'k', lw=3)
    plt.plot(wn, pure_spectra.T, 'b')
 
    plt.xlabel("Wavenumber (cm$^{-1}$)", fontsize=14)
    plt.ylabel("Raman spectra", fontsize=14)
 
    plt.show()

Исходные чистые спектры показаны черным цветом, а спектры, восстановленные с помощью подхода MCR, - синим. Кроме того, в массиве mcrar.C_ доступны реконструированные концентрации, которые можно сравнивать с реальными концентрациями.

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

Спасибо за прочтение и до следующего раза,
Дэниел.

Разрешение многомерной кривой: введение

Спектры синтетических смесей

Нахождение чистых спектров по концентрациям

Нахождение концентраций по чистым спектрам

Разрешение многомерной кривой

Рекомендации

Вопросы по теме