Об использовании Visual Studio Code для учебных примеров

Visual Studio Code — редактор на основе платформы Electron, который можно использовать как минималистичную IDE. Мне требовалась поддержка C++17 и переносимость всего «на флешке». Что меня приятно удивило, так это «шустрость» VSCode.

Итак, порядок действий, предпринятых мной был таков:

  1. Скачать VSCode в архиве (в моём случае это была версия 1.27.1 для Windows/x64). Распаковать его, допустим, в папку vscode. Внимание! Путь к vscode не должен содержать не-ASCII символы и пробелы.
  2. Создать пустую папку vscode/data. Там будут плагины и общие настройки. При обновлении VSCode можно просто распаковать новую версию и переместить туда свою папку data из старой.
  3. Распаковать сборку MinGW рядом с vscode, допустим, в папку mingw (mingw/bin содержит g++.exe). В качестве таковой я брал сборку отсюда, там уже есть Boost и ряд других библиотек. Я к ним ещё добавляю Eigen и кое-что своё. Добавляю, просто складывая соответствующие папки в mingw/include.
  4. Создать рядом с vscode пустую папку для примеров. Допустим, cpp.
  5. Запустить VSCode (vscode/code.exe) и установить плагин C/C++ for Visual Studio Code от Microsoft (с Microsoft Intellisense).
  6. Открыть в VSCode папку cpp.
  7. Настройки создаются в виде JSON-файлов в папке cpp/.vscode. Ниже я опишу возможное их содержимое. Их можно создать непосредственно в VSCode, заполнить в нём же, и затем перезапустить VSCode.
  8. После того, как файлы настроек созданы, и VSCode перезапущен, следует протестировать его. Создать cpp/test.cpp с произвольным содержимым на C++ (например, helloworld). Terminal/Run build task… или Ctrl+Shift+B должно запускать сборку и завершаться без ошибок. F5 запускает сборку и результат в отладчике (gdb). Заголовочные файлы должны успешно открываться (правой кнопкой мыши и Go to definition). Intellisense должен работать адекватно.
  9. Полученный «пакет» в составе vscode, mingw, cpp можно упаковать в архив и развернуть на другой системе или просто скопировать на внешний носитель и запускать с него.

Итак, файлы настроек.

Плагин C/C++ использует файл c_cpp_properties.json.
Его содержимое в моём случае:

{
    "configurations": [
        {
            "name": "MinGW",
            "intelliSenseMode": "gcc-x64",
            "compilerPath": "${workspaceFolder}\\..\\mingw\\bin\\g++.exe",
            "cppStandard": "c++17",
            "cStandard": "c11",
            "includePath": [
                "${workspaceFolder}",
                "${workspaceFolder}\\..\\mingw\\include"
            ],
            "defines": [
                "__cplusplus=201703L"
            ]
        }
    ],
    "version": 4
}

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

Для сборки требуется создать «задачу сборки». В простейшем случае это компиляция одной единицы трансляции (cpp-файла) сразу в исполняемый файл (в моём случае всегда один и тот же с названием lab.exe).
Описание задач сборки располагается в файле tasks.json. Его содержимое в моём случае:

{
    // See https://go.microsoft.com/fwlink/?LinkId=733558
    // for the documentation about the tasks.json format
    "version": "2.0.0",
    "tasks": [
        {
            "label": "build lab",
            "type": "shell",
            "command": "${workspaceFolder}\\..\\mingw\\bin\\g++",
            "args": [
                "-g", "-Og", // "debug" build
                //"-O3",
                "${file}", 
                "-march=native", "-Wall", "-pedantic",
                "-std=c++17",
                "-olab.exe"
            ],
            "group": {
                "kind": "build",
                "isDefault": true
            }
        }
    ]
}

Наконец, описание запуска отладчика по F5 располагается в launch.json.
Его содержимое в моём случае:

{
    // Use IntelliSense to learn about possible attributes.
    // Hover to view descriptions of existing attributes.
    // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387
    "version": "0.2.0",
    "configurations": [
        {
            "name": "(gdb) Launch",
            "type": "cppdbg",
            "request": "launch",
            "program": "${workspaceFolder}\\lab.exe",
            "args": [],
            "stopAtEntry": false,
            "cwd": "${workspaceFolder}",
            "environment": [],
            "externalConsole": true,
            "MIMode": "gdb",
            "miDebuggerPath": "${workspaceFolder}\\..\\mingw\\bin\\gdb.exe",
            "setupCommands": [
                {
                    "description": "Enable pretty-printing for gdb",
                    "text": "-enable-pretty-printing",
                    "ignoreFailures": true
                },
                {
                    "description": "Enable break on exception for gdb",
                    "text": "catch throw",
                    "ignoreFailures": true              
                }
            ],
            "preLaunchTask": "build lab"
        }
    ]
}

Вторая команда catch throw активирует перехват исключений отладчиком gdb.
Вообще, к сожалению, я не могу сказать, что с отладкой всё хорошо: есть ряд странностей. Например, точка останова, поставленная в одном месте, может сработать ниже (как будто поставленная в другом месте — может быть, это можно исправить, если заменить -Og на -O0). Точку останова нельзя поставить внутри шаблона (не срабатывает вообще). Не даёт ставить новые точки останова во время отладки. Если во время отладки нажать паузу, то продолжение может не срабатывать (останется только остановить исполнение).

Для меня проблемы с отладчиком не являются критичными, так как я привык к отладке через проверку условий и протоколирование. Но тем, кто привык к отладчику в «нормальной» Visual Studio такое, конечно, будет не по нраву.

UPD. Если всё это запускать на Windows 10, то в качестве терминала может быть по умолчанию выбран PowerShell. Почему-то он «не понимает» относительные пути через .. и не может запустить компилятор. В этом случае рекомендуется заменить его на cmd.exe (в VSCode можно выбрать терминал).

UPD2. Замена флага оптимизации -Og на -O0 всё же помогает отладчику, по крайней мере, точки останова не «убегают» с тех мест, куда они были поставлены.

Реклама

Граница Парето в 2D: STL-версия

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

С тех пор я использовал данную задачу как задачу по программированию на первом курсе (на алгоритмы STL). Всё дело в том, что этот алгоритм очень просто записывается с помощью STL. Требуется применить два алгоритма (sort и unique) и адаптер компаратора (not_fn или самописный, например, переставляющий аргументы местами). Бывают простенькие задачи, к которым невозможно адекватно применить STL без написания своих компонент, а бывает и наоборот. Итак, код.

#include <algorithm>
#include <functional>

template 
  <
    class RanIt, 
    class Comp1, 
    class Comp2
  >
RanIt remove_dominated
  (
    RanIt from,
    RanIt to,
    Comp1 compare1,
    Comp2 compare2
  )
{
  std::sort(from, to, std::not_fn(compare1));
  return std::unique(from, to, std::not_fn(compare2));
}

////////////////////////////
#include <string>
#include <vector>
#include <iostream>
using namespace std;

bool test_remove_dominated()
{
  vector<string> input
  {
    "look",
    "after",
    "the",
    "raven",
    "star"
  };
  
  input.erase(
    remove_dominated(
      input.begin(), input.end(),
      std::less<>{},
      [](auto & a, auto & b)
      {
        return a.size() < b.size();
      }),
    input.end()
    );
    
  sort(input.begin(), input.end());
  return input == vector<string>
    {
      "raven",
      "star",
      "the"
    };
}

int main()
{
  cout << test_remove_dominated() << '\n';
  return 0;
}

В примере два критерия сравнения строк: лексикографический и по длине. Так «raven» доминирует «look» (и лексикографически, и по длине), но не «star».

О задаче разложения числа в сумму

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

Дано: x \in \mathbb R, n \in \mathbb N, \{\underline{x}_i\}_{i=1}^n \subset \mathbb R, \{\overline{x}_i\}_{i=1}^n \subset \mathbb R, \forall i \in \overline{1,n}\, (\underline{x}_i \leqslant \overline{x}_i).

Найти \{x_i\}_{i=1}^n такие, что \sum_{i=1}^n{x_i} = x и \forall i \in \overline{1,n}\, (\underline{x}_i \leqslant x_i \leqslant \overline{x}_i), или определить, что таких x_i не существует.

Задача легко решается, если догадаться сделать замену переменных:

y_i = x_i - \underline{x}_i, \quad \overline{y}_i = \overline{x}_i - \underline{x}_i.

y = \sum\limits_{i=1}^n{y_i} = \sum\limits_{i=1}^n{(x_i - \underline{x}_i)} = x - \sum\limits_{i=1}^n{\underline{x}_i}.

Она позволяет перейти к более простой постановке:

Дано: y \in \mathbb R, n \in \mathbb N, \{\overline{y}_i\}_{i=1}^n \subset \mathbb R, \forall i \in \overline{1,n}\, (0 \leqslant \overline{y}_i).

Найти \{y_i\}_{i=1}^n такие, что \sum_{i=1}^n{y_i} = y и \forall i \in \overline{1,n}\, (0 \leqslant y_i \leqslant \overline{y}_i), или определить, что таких y_i не существует.

Сразу видно, что решение существует тогда и только тогда, когда 0 \leqslant y \leqslant \sum_{i=1}^n{\overline{y}_i} и получить его можно простым линейным жадным алгоритмом, вычитая из y на каждом шаге максимально возможный очередной недостаток y_i = \min\{\,y - \sum_{j=1}^{i-1}{y_j}, \overline{y}_i\,\} (естественно, y_1 = \min\{\,y, \overline{y}_1\}).
Если мы быстро исчерпаем y, то последние y_i могут быть просто нулями.

Приведу код решения данной задачи на C++ (аж с ноября не публиковал здесь кода C++!).

#include <cassert>
#include <iterator>
#include <algorithm>
/**
 * @brief Decompose a number into a sum of non-negative terms.
 * @param value the value to be decomposed
 * @param ubound non-negative upper bounds for decomposition terms
 * @param decomp decomposition terms, of the same size as ubound!
 * @return success (whether a decomposition exists)
 */
template 
  <
    class NT,       // number type
    class NTArrayU, // number type array for ubound
    class NTArrayD  // number type array for decomp
  >
bool scalar_sum_decompose
  (
    NT               value,
    NTArrayU const & ubound,
    NTArrayD       & decomp
  )
{
  if (value < NT{})
    return false;
  
  using std::size;
  auto const n = size(ubound);
  assert(n == size(decomp));
  
  NT s {}; // accumulated sum.
  for (auto i = n-n; i < n; ++i)
  {
    assert(!(ubound[i] < NT{}));
    NT const term = std::min<NT>(value - s, ubound[i]);
    decomp[i] = term;
    value -= term;
  }
  
  return value == NT{};
}

(Код «auto i = n-n» есть маленький костыль ввиду нежелания городить извлечение типа размера.)
Теперь несложно написать и обёртку, решающую первую задачу.

#include <vector>
#include <numeric>
#include <functional>
/**
 * @brief Decompose a number into a sum of bounded terms.
 * @param value the value to be decomposed
 * @param lbound lower bounds for decomposition terms
 * @param ubound upper bounds for decomposition terms
 * @param decomp decomposition terms, it must have the same size as ubound and lbound!
 * @return success (whether a decomposition exists)
 */
template 
  <
    class NT,       // number type
    class NTArrayL, // number type array for lbound
    class NTArrayU, // number type array for ubound
    class NTArrayD  // number type array for decomp
  >
bool scalar_sum_decompose
  (
    NT               value,
    NTArrayL const & lbound,
    NTArrayU const & ubound,
    NTArrayD       & decomp
  )
{
  using std::size;
  auto const n = size(ubound);
  assert(size(lbound) == n && size(decomp) == n);
  
  using std::begin;
  using std::end;
  std::transform(
    begin(ubound), end(ubound), begin(lbound),
    begin(decomp), std::minus<>{});
  NT const lsum = std::accumulate(
    begin(lbound), end(lbound), NT{});
  std::vector<NT> buffer(n);
  if (!scalar_sum_decompose(value - lsum, decomp, buffer))
    return false;
  std::transform(
    begin(buffer), end(buffer), begin(lbound),
    begin(decomp), std::plus<>{});
  return true;
}

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

Имя/enmen

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

Начнём с предлога в : в < въ < *vъn/*vьn (-n сохранено в некоторых производных: внутрь = вн-утрь (ср. утроба), вниди = вн-иди; протетическая в- как в вострый, вогнь, восемь) < *en < *h₁én (см. там также другие производные).

Сюда, может быть, также вон (вънъ), вне (вънѣ) (у Фасмера по-другому).
Последняя ассоциация у меня возникла через другую ассоциацию с местоимением он < *h₁ónos («тот, удалённый от меня и от тебя»; «некий», «какой-то»). Ещё варианты этого местоимения по ссылке: *h₁énos (*h₁én-o-s?), *h₁nós.

О, эта последняя форма прекрасна. Она мне сразу напомнила слово *h₁nómn̥ «имя» (см. также спор в обсуждении о возможном корне этого слова).

Обычно предполагается, что *-mn̥ (-men) — суффикс отглагольного существительного, тогда *h₁nó-mn̥ должно быть производным некоего глагола со значением «давать имя» или «звать», ср. прозвище от прозвать. Думаю, можно найти множество подобных примеров в существующих языках. Попытка истолковать данное слово на индоевропейской почве отягощается сравнением с материалом других языковых семей, например, с прауральским *nime «имя».

Балтослав. реконструкция *inˀmen может быть объяснена влиянием форм косвенных падежей (форм вроде *h₁n̥mén, см. таблицу склонения по ссылке выше). Действительно, примеры похожих замен имеются: например, слова язык и земля (впрочем, это, по сути, упрощения исходных сложных слов).

Ввиду древнепрусского emmens ЭССЯ (т.8, с.228) предлагает понимать «имя» как *en-men «влагаемое» > «имя» и считает его параллельным образованию *ano-men «возлагаемое» > «имя». Наличие различных (но похожих) параллельных образований объясняет невозможность единой непротиворечивой реконструкции. Уральские формы там объяснены как ИЕ заимствование (это утверждение не будет таким уж сильным, если принять возможную датировку позднего прауральского языка 4000 лет назад).

Можно потеоретизировать, что форма *h₁nós даёт возможность некоего подобного толкования, но уже на основе указательного местоимения. Вообще, наверно, не будет ошибкой считать указательные местоимения наиболее примитивной частью речи. И как таковые, они могли служить на некоем древнем этапе базисом для других слов, например, для предлогов направления («в» и «на»).

Наконец, интересно сходство слав. имя с глаголом иметь (атематический ять < *imtei < *h₁em-). Это можно было бы назвать случайностью, но последний корень сравнивают с корнем *nem-, имеющим аналогичное значение. Отглагольное производное с корнем *nem/nom могло бы значить «данное (имя)».

Чё и что

В современном русском языке слово чё считается просторечным вариантом что. Тем не менее, оно является прямым продолжением индоевропейского вопросительного слова. Для пояснения приведу таблицу изменений по этапам (включая также местоимения 3л.).

Этап что кто он, оно, она
ПИЕ *kʷid *kʷos *éy, *Hyós, *Hyéh₂
сатемизация *kid *kos *ey, *yos, *ya
славизация *kь *kъ *jь, *je, *ja
первая палатализация *чь *къ (ѥго), (ѥго), (ѥѩ)
падение еров че **ко (им.п. заменён)

Последний исторический этап: переход подударного е в ё в русском языке (к XVIII веку?).

Видимо, из-за того, что формы именительного падежа стали в процессе изменения слишком короткими и омонимичными (например, къ — предлог, и — союз), они были заменены, соответственно, на что (чьто, чь+то), кто (къто, чъ+то), он, оно, она.

Другие составные местоимения: чей (чь+и), кой (къ+и), где (къ+дѣ/де), когда (возможно, *кого года, т.е. в какое время) и т.д. Комбинирование местоимений в славянском дало комбинаторный взрыв числа их вариантов, из которых в современных языках закрепились не все (порой, разные слова для одинаковых значений).

О вычислении матричной экспоненты

1. Экспонента квадратной матрицы \mathbf A формально определяется как сумма ряда:

e^{\mathbf A} \triangleq \mathbf I + \sum\limits_{i=1}^{\infty} \dfrac{\mathbf A^i}{i!},

где \mathbf I есть единичная матрица соответствующего размера.

Из этого определения несложно получить ряд свойств экспоненты, некоторые из которых перечислены ниже:

  1. Экспонента нулевой матрицы есть единичная матрица.
  2. Транспонирование (а также эрмитово сопряжение) перестановочно с экспоненцированием: e^{\mathbf A^\top} = (e^{\mathbf A})^\top.
  3. Если \mathbf A симметрична (эрмитова), то и e^{\mathbf A} симметрична (эрмитова).
  4. Если \mathbf {A B} = \mathbf {B A}, то e^{\mathbf A} \mathbf B = \mathbf B e^{\mathbf A} и e^{\mathbf A} e^{\mathbf B} = e^{\mathbf B} e^{\mathbf A} = e^{\mathbf A + \mathbf B}. В частности, e^{a \mathbf A} e^{b \mathbf A} = e^{(a + b)\mathbf A}.
  5. Матрица e^{\mathbf A} всегда обратима: e^{\mathbf A} e^{-\mathbf A} = \mathbf I.
  6. Для любой обратимой матрицы \mathbf P верно e^{\mathbf{P A P}^{-1}} = \mathbf P e^{\mathbf A} \mathbf P^{-1}.
  7. Экспонента диагональной матрицы \mathbf A = \mathrm{diag}(a_{ii}) также является диагональной матрицей e^{\mathbf A} = \mathrm{diag}(e^{a_{ii}}).
  8. Последнее свойство обобщается на блочно-диагональные матрицы. Пусть \mathbf A = \mathrm{diag}(\mathbf A_k), где каждая \mathbf A_k есть квадратная матрица. Тогда экспонента также является блочно-диагональной матрицей: e^{\mathbf A} = \mathrm{diag}(e^{\mathbf A_k}).

Таким образом, если матрица \mathbf A диагонализируема, т.е. если найдётся такая матрица \mathbf P, что \mathbf A = \mathbf {P D P}^{-1}, где \mathbf D = \mathrm{diag}(\lambda_i) — диагональная матрица (матрица собственных значений), то после выполнения диагонализации (eigendecomposition) \mathbf A вычислить её экспоненту не составит труда.

К сожалению, не все матрицы диагонализируемы, хотя для матриц над полем комплексных чисел множество недиагонализируемых матриц есть множество меры нуль. Поэтому, кстати, даже экспоненту вещественнозначной матрицы может быть удобнее вычислять в комплексных числах. Те, кто хорошо учил линейную алгебру, вспомнят о «жордановой форме», которая существует для любой матрицы и представляет собой как раз блочно-диагональную матрицу, составленную из блоков, называемых «жордановы клетки». Это означает, что если у нас есть представление \mathbf A = \mathbf{P J P}^{-1}, где \mathbf J = \mathrm{diag}(\mathbf J_k), и если мы умеем вычислять экспоненту жордановой клетки, то мы опять же без проблем вычислим экспоненту исходной матрицы. Ниже я коснусь вопроса вычисления экспоненты жордановой клетки.

Особняком стоит свойство, известное как «формула Якоби». Его можно применять, например, для проверки качества вычисления матричной экспоненты:

\det e^{\mathbf A} = e^{\mathrm{tr}\,\mathbf A}.

2. Вычисление экспоненты прямым суммированием ряда традиционно считается «плохим методом». На самом деле, главный его недостаток заключается в низкой скорости вычислений, но не в точности. Важно организовать этот процесс правильно! В своей практике я применял следующие приёмы:

  • Компенсационное суммирование. Например, алгоритм Ноймайера.
  • Вместо вычисления очередной степени \mathbf A^i домножением, я вычислял степень каждый раз «заново» (до 2\log_2{i} + 1 матричных умножений).
  • Единичная матрица добавляется в конце (опционально: не нужно её добавлять, если вычисляется e^{\mathbf A} - \mathbf I ).
  • Суммирование продолжается до тех пор, пока очередное слагаемое не обнулится из-за слишком большого делителя-факториала (на двойной точности обычно достаточно несколько десятков слагаемых).

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

3. Экспонента жордановой клетки. Жорданова клетка \mathbf J (в случае комплексных чисел) есть квадратная матрица, на главной диагонали которой стоит некое число \lambda (любое значение, но одно во всех позициях диагонали), на первой наддиагонали стоят единицы, а все прочие элементы — нули:

\mathbf J = \begin{bmatrix}  \lambda & 1 & 0 & \ldots & 0 \\  0 & \lambda & 1 & \ldots & 0 \\  & & \cdots & & \\  0 & 0 & \ldots & \lambda & 1 \\  0 & 0 & \ldots & 0 & \lambda  \end{bmatrix}.

Заметим, что \mathbf J = \mathbf \Lambda + \mathbf N, где \mathbf \Lambda — шаровая матрица (диагональная с одинаковыми значениями на диагонали), а \mathbf N — нильпотентна, поскольку \mathbf N^n = \mathbf 0 (здесь и ниже n — размер \mathbf J). Более того, легко заметить, что домножение \mathbf N^i на \mathbf N попросту сдвигает единичную наддиагональ на одну позицию вверх. Поэтому имеем

e^{\mathbf N} = \mathbf I + \sum\limits_{i=1}^{n-1}\dfrac{\mathbf N^i}{i!} = \begin{bmatrix}  1 & 1 & 1/2 & 1/6 & \ldots & 1/(n - 1)! \\  0 & 1 & 1 & 1/2 & \ldots & 1/(n-2)! \\  & & & \cdots & & \\  0 & 0 & \ldots & 0 & 1 & 1 \\  0 & 0 & \ldots & 0 & 0 & 1  \end{bmatrix}.

Ну и по свойствам экспоненты и шаровых матриц имеем

e^{\mathbf J} = e^{\mathbf \Lambda + \mathbf N} = e^{\mathbf \Lambda} e^{\mathbf N} = e^\lambda e^{\mathbf N}.

В случае наличия скалярного множителя t:

e^{\mathbf J t} = e^{\mathbf \Lambda t + \mathbf N t} = e^{\lambda t} e^{\mathbf N t},

где

e^{\mathbf N t} = \begin{bmatrix}  1 & t & t^2/2 & t^3/6 & \ldots & t^{n-1}/(n - 1)! \\  0 & 1 & t & t^2/2 & \ldots & t^{n-2}/(n-2)! \\  & & & \cdots & & \\  0 & 0 & \ldots & 0 & 1 & t \\  0 & 0 & \ldots & 0 & 0 & 1  \end{bmatrix}.

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

4. О вычислении матричной экспоненты (и иных функций, определяемых как сумма ряда) с помощью комбинирования разложения Шура, рекуррентных соотношений Сильвестра и суммировании ряда см. Davies, Higham. A Schur-Parlett Algorithm for Computing Matrix Functions.

Об аналитическом решении линейных ОДУ

1. Пусть дана однородная система (задача Коши):

\dot{\mathbf x} = \mathbf A \mathbf x, \quad \mathbf x(t_0) = \mathbf x_0,

где \mathbf x(t) \in \mathbb R^n — фазовый вектор системы, \mathbf A — постоянная n\times n-матрица над полем \mathbb R. Цель — найти функцию \mathbf x(t), t \geqslant t_0, удовлетворяющую указанной системе. Запись \dot{\mathbf x} обозначает производную \mathbf x'(t).

Как и в случае n = 1 решение можно получить через экспоненту.
В нашем случае — матричную экспоненту, определяемую формально как сумму ряда. Вычисление матричной экспоненты, особенно для больших n, есть отдельная большая тема, которой в данной статье я не намерен касаться подробно. Ниже приведены формулы для прямого вычисления матричной экспоненты в случае n = 2.

Итак, сразу решение:

\mathbf x(t) = e^{\mathbf A(t - t_0)} \mathbf x_0.

Заметим, что, как и в скалярном случае, мы можем продифференцировать функцию e^{\mathbf A t} по t как матрицу покомпонентно, получив:

(e^{\mathbf A t})' = e^{\mathbf A t} \mathbf A = \mathbf A e^{\mathbf A t}.

\mathbf A, \mathbf A^{-1} и e^{\mathbf A t} являются перестановочными матрицами, что легко показать, раскрыв экспоненту в сумму ряда.

2. Пусть теперь дана неоднородная система:

\dot{\mathbf x} = \mathbf A \mathbf x + \mathbf f, \quad \mathbf x(t_0) = \mathbf x_0,

где \mathbf f \colon \mathbb R \to \mathbb R^n. Выполнив замену

\mathbf z(t) = e^{\mathbf A(t_0 - t)} \mathbf x(t),

исключим фазовую переменную из правой части уравнения:

\dot{\mathbf z}(t) = e^{\mathbf A(t_0 - t)} \mathbf f(t), \quad \mathbf z(t_0) = \mathbf x_0,

откуда легко получить окончательную формулу для решения исходной неоднородной системы:

\mathbf x(t) = e^{\mathbf A(t - t_0)} \left( \mathbf x_0 + \int_{t_0}^t {e^{\mathbf A(t_0 - \tau)} \mathbf f(\tau)} \,d\tau\right).

3. Пусть \mathbf f(t) \equiv \mathbf u = const. Пусть также \det\mathbf A \neq 0. Тогда имеем первообразную:

(e^{\mathbf A t} \mathbf A^{-1})' = e^{\mathbf A t}.

После подстановки первообразной в интеграл формула решения из предыдущего пункта обретает вид:

\mathbf x(t) = e^{\mathbf A(t - t_0)} (\mathbf x_0 + \mathbf A^{-1} \mathbf u) - \mathbf A^{-1} \mathbf u.

4. Предположим теперь, что матрица \mathbf A является функцией времени. Оставим также функцию \mathbf f. Можно воспользоваться тем же приёмом, что и в п.2, введя замену:

\mathbf z(t) = e^{\int_t^{t_0} {\mathbf A(\tau)}\,d\tau} \mathbf x(t).

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

\forall t_1, t_2 \in [t_0, t] \left(\,\mathbf A(t_1) \mathbf A(t_2) = \mathbf A(t_2) \mathbf A(t_1)\,\right).

Тогда можно записать решение системы в форме:

\mathbf x(t) = e^{\int_{t_0}^t {\mathbf A(\tau)}\,d\tau} \left( \mathbf x_0 + \int_{t_0}^t {e^{\int^{t_0}_\tau {\mathbf A(s)}\,ds} \mathbf f(\tau)} d\tau\right).

5. Формулы для матричной экспоненты размерности 2. Пусть дана матрица

\mathbf A = \begin{bmatrix} a & b \\ c & d  \end{bmatrix}

Положим D = 4 bc + (a - d)^2. Имеем теперь три случая.

Случай D > 0.
Положим q = \sqrt{D}, C = \cosh\frac{q}{2}, S = \sinh\frac{q}{2}. Тогда

e^{\mathbf A} = \frac{1}{q} e^\frac{a+d}{2}  \begin{bmatrix} qC + (a - d) S & 2bS \\                             2cS & qC - (a - d)S  \end{bmatrix}

Случай D < 0.
Формула та же, что и в случае D > 0, если положить q = \sqrt{-D}, C = \cos\frac{q}{2}, S = \sin\frac{q}{2}.

Случай D = 0.

e^{\mathbf A} = e^\frac{a+d}{2}  \begin{bmatrix} 1 + \frac{a - d}{2} & b \\                             c & 1 - \frac{a - d}{2}  \end{bmatrix}

Для множества частных случаев можно вывести более простые формулы. Например, для диагональной матрицы b = c = 0:

e^{\mathbf A} = \begin{bmatrix}  e^a & 0 \\ 0 & e^d \end{bmatrix}

Для верхнетреугольной матрицы c = 0:

e^{\mathbf A} = \begin{bmatrix}  e^a & b\frac{e^a - e^d}{a - d}   \\ 0 & e^d \end{bmatrix}

О методе главных компонент

На днях пришлось столкнуться с «методом главных компонент» (principal component analysis, PCA). Ниже даю обретённое мною понимание оного.

Дано: набор точек { pi ∈ ℝd | i∈1:n }, dn.

Предположение: точки получены аффинным преобразованием реализации случайного процесса y: pi = μ + yiW, где μ ∈ ℝd, W — некая матрица, а yi удовлетворяют многомерному нормальному распределению с нулевым матожиданием и единичной матрицей ковариации.

Задача: найти μ, W и yi.

Замечание 1. Размерность yi не обязана совпадать с размерностью pi.

Замечание 2. В качестве оценки матожидания μ возьмём арифметическое среднее pi: μ = Σi pi / n. Обозначим: xi = pi − μ.

Замечание 3. Так как покоординатно независимое нормальное распределение порождает «шарообразное» скопление точек yi, то следует ожидать, что xi формируют скопление, визуально напоминающее повёрнутый d-осный эллипсоид. Иными словами, можно искать матрицу W как комбинацию анизотропного масштабирования и поворота: W = TR, где T может иметь ненулевые элементы только на главной диагонали, а R ортогональна.

Через X обозначим n×d-матрицу, составленную из строк xi.

«Поточечная» матрица ковариации XX имеет размеры n×n. «Покоординатная» матрица ковариации XX имеет размеры d×d. Обе эти матрицы симметричны, но в общем случае не диагональны.

Через Y обозначим n×r-матрицу, составленную из строк yi. По предположению имеем YY = 1n или YY = 1r. Выберем пока первый вариант. Из него немедленно следует rn. Логично положить r = n, тогда Y ортогональна.

Итак, X = YTR, где Y и R ортогональны, а T имеет ненулевые элементы только на главной диагонали. Данное разложение матрицы X есть ни что иное, как разложение сингулярных значений (singular value decomposition, SVD).

Возьмём любую реализацию SVD и построим разложение: X = USV*. В общем случае матрицы U и V унитарны, звёздочкой обозначено эрмитово сопряжение. В нашем случае они ортогональны и сопряжение есть транспонирование, поэтому можно положить:

Y = U, T = S, R = V, W = SV.

Замечание 4. На диагонали S стоят сингулярные значения, в нашей задаче имеющие геометрический смысл длин полуосей эллипсоида. Они неотрицательны и по построению отсортированы по величине от больших к меньшим. В конце могут стоять нули, означающие, что по соответствующим направлениям (которым отвечают столбцы V) эллипсоид вырожден, и эти измерения можно отбросить. С практической точки зрения может быть разумным отбросить все измерения, для которых полуось (сингулярное значение) меньше некоторой заданной величины. WW = SS — диагональная матрица, составленная из квадратов сингулярных значений.

Шаг 1. A = XV. Так как V есть матрица, обратная к R, то A = YT. Это эллипсоид, повёрнутый так, чтобы его оси совпали с осями координат стандартного базиса, причём разброс данных на каждой следующей оси не больше, чем на предыдущей. Кому-то этого результата уже может быть достаточно.

Шаг 2. Пусть задано ρ ≥ 0 — предельная «чувствительность» по измерению. Исходная полная матрица S может быть неквадратной и содержать нулевой блок под квадратом d×d, здесь я его проигнорирую и положу S = diag(sj), 1 ≤ jd.

Пусть r = max{ j | sj > ρ }. Положим B = A[1:r], т.е. первые r столбцов матрицы A. Строки bi матрицы B есть точки в пространстве уменьшенной размерности r ≤ d.

Шаг 3. Составим усечённую псевдообратную матрицу S+ = diag(1/sj), 1 ≤ jr.
Положим Z = BS+. Так как S+ есть матрица, обратная к усечённой T, то мы получим, что Z = Y[1:r], т.е. первые r столбцов Y. Если мы уже вычислили U, то можно просто взять из неё первые r столбцов: Z = U[1:r]. Имеем ZZ = 1r.

Точки zi имеют уменьшенную размерность r ≤ d и отмасштабированы по осям координат.

Гладкая интерполяция между синусоидой и прямоугольной волной

Представим, что есть задача получить нечто среднее между синусоидой

\mathrm{Sin}\,(t) \triangleq \sin{2\pi t}

и прямоугольной волной

\mathrm{Rect}\,(t) \triangleq \mathrm{sgn}\,\mathrm{Sin}\,(t).

Очевидно, что банальная поточечная линейная интерполяция с параметром \alpha \in [0, 1]:

(1-\alpha)\,\mathrm{Sin}\,t + \alpha\,\mathrm{Rect}\,t

не является функцией, всюду дифференцируемой по t. Впрочем, оказывается, легко получить гладкую «промежуточную» функцию, использовав арктангенс:

R_\alpha\,(t) \triangleq \dfrac{\mathrm{arctg}\,(\alpha\,\mathrm{Sin}\,(t))}{\mathrm{arctg}\,\alpha}.

Здесь уже \alpha \in (0, +\infty).
Легко видеть, что

\lim\limits_{\alpha\to 0} R_\alpha\,(t) = \mathrm{Sin}\,(t),\quad \lim\limits_{\alpha\to +\infty} R_\alpha\,(t) = \mathrm{Rect}\,(t).

Графики R_\alpha\,(t), t \in [0, 2], \alpha \in \{\,1, 10, 1000\,\}:

Один алгоритм раскраски графа II

Недавно я переписал алгоритм субоптимальной правильной вершинной раскраски графа на Python. Решил выложить его сюда.

Вершина графа отождествляется с её индексом. Граф задаётся индексируемой последовательностью (tuple или list) наборов соседей (массив списков соседей). Раскраска — массив (list в терминах Python) цветов.

Итак, функция раскраски у меня имеет следующий вид:

## Построить правильную вершинную раскраску графа.
## \param graph -- граф в виде массива множеств номеров смежных вершин
## \param forbidden -- массив множеств запрещённых цветов для каждой вершины
## \param color -- массив цветов вершин
## \return массив (новых) цветов вершин
def coloring(graph, forbidden=None, color=None):
    if forbidden is None: forbidden = []
    if color is None: color = []
    # Предусловия.
    assert len(forbidden) == 0 or len(forbidden) == len(graph)
    assert len(color) == 0 or len(color) == len(graph)

    if len(forbidden) == 0:
        forbidden = [frozenset()]*len(graph)
    #...

Передав раскраску color, можно задать желаемые цвета, а алгоритм разрешит конфликты, если они есть. Если color уже содержит правильную раскраску, то она и будет возвращена функцией coloring. Массив forbidden позволяет запретить раскрашивать какие-то вершины в какие-то цвета: forbidden[v] есть множество запрещённых цветов (числовых меток) для вершины с индексом v.

Всё, что описано ниже, определено внутри функции coloring.

Если пользователь не передал некую начальную раскраску color, то построим её в два цвета поиском в глубину.

    # Раскраска в два цвета поиском в глубину.
    def coloring2():
        color2 = [0]*len(graph)
        def dfsColoring2(u, c=1):
            if color2[u] == 0:
                for d in range(c, len(graph)+1):
                    if d not in forbidden[u]:
                        color2[u] = d
                        break
                c = 3 - c
                for v in graph[u]:
                    dfsColoring2(v, c)
    
        for start in range(0, len(graph)):
            dfsColoring2(start)
        return color2

    # Построить начальную раскраску, если требуется.
    if len(color) == 0:
        color = coloring2()

Теперь собственно алгоритм. Он заключается в повторении пары действий: выявлении конфликтов цветов (listConflicts) и разрешении конфликтов цветов (resolveConflicts), пока все конфликты не будут разрешены.

    # Перечислить текущие конфликты цветов.
    def listConflicts():
        conflicts = []
        for u in range(0, len(graph)):
            nc = [color[v] for v in graph[u]]
            conflict_rank = nc.count(color[u])
            if conflict_rank > 0:
                nc = frozenset(nc) | forbidden[u]
                new_color = max(nc) + 1
                for c in range(1, new_color-1):
                    if c not in nc:
                        new_color = c
                        break
                conflict = (conflict_rank, new_color, u)
                conflicts.append(conflict)
        conflicts.sort(
            key=lambda c: (c[0], -c[1], c[2]), reverse=True)
        return conflicts
    
    # Разрешить конфликты цветов (потенциально не все).
    def resolveConflicts(conflicts):
        closed = set()
        for _, new_color, u in conflicts:
            if u in closed:
                continue
            color[u] = new_color
            closed |= set(graph[u])
    
    # Собственно алгоритм: начиная с некоторой раскраски,
    # повторять обнаружение и разрешение конфликтов, 
    # пока список конфликтов не опустеет.
    while True:
        conflicts = listConflicts()
        if len(conflicts) == 0:
            return color
        resolveConflicts(conflicts)

Цвет 0 не используется.

Пример:

g = ((1,2), (0,2), (0,1,3,6), (2,4,5), (3,5,6), (3,4), (2,4))
c = coloring(g)
print(c)
> [1, 2, 3, 2, 1, 3, 2]

Раскрашенный граф: