Добрый день, группа 312!

План занятия

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

Материалы к занятиям: https://maxxk.github.io/programming-semester-5/
email:

В дисплейных классах рекомендуется просматривать в браузере Firefox.
В нём установлено расширение NoScript, обратите внимание на инструкцию, иначе значительная часть сайтов не будет работать.

Числа с плавающей точкой

Плавающая точка — метод представления подмножества рациональных чисел (бывает тип комплексных чисел с плавающей точкой):

$
x = (-1)^{{\color{#0CC}{sign}}} · 1.\overline{{\color{#F00} {significand}}} · 2^{{\color{#0C2}{exponent - bias}}}$

32-битное число (float):
floating-point

Смещенный порядок (biased exponent): p = exp − 127

Нормированная мантисса (normalized significand): ведущий бит опускается и считается равным единице

На картинке мантисса — $ \color{#444}{1}.01000…0 $

Подробнее
Подробнее на Wikipedia (и далее — про стандарт IEEE 754)

Особые значения

Мантисса — двоичная дробь. Не всякая десятичная дробь представляется как конечная двоичная.
1/5 (0.2) записывается в double как:
$0.2000000000000000\color{red}{11102230246251565404236316680908203125}$
⇒ для финансовых приложений нужно использовать десятичную запись с фиксированной точкой и строго заданными правилами округления
Денормализованные числа — порядок = 0, но мантисса не может быть записана в нормированном виде — число записывается в денормализованном виде, что вызывает потерю точности значащих цифр.
NaN (не-число) (все биты поряка установлены) — возвращается в операциях, результат которых не определён (0/0, $\sqrt{-1}$)
Infinity (±∞) — возвращается при делении на 0 или переполнении (в дисплейных классах вместо ∞ — исключение)

Исключения

Floating Point Exception (SIGFPE):

некорректная операция ($\sqrt{-1}$)
деление на 0 (частая причина — n = 1 в x/(n − 1))
переполнение (частая причина — бесконечный цикл/невыполнение условия выхода из цикла, например, сравнение чисел с плавающей точкой оператором ==)
underflow (частая причина — обращение к неинициализированному элементу массива или выход за границы массива)

Segmentation Fault (SIGSEGV):

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

Постановка задачи

𝐌_n — кольцо матриц n×n над ℂ (со сложением и умножением на константу можно рассматривать как векторное пространство над ℂ)
I — единичная матрица (E — матрица погрешности)

Решение линейной системы:
найти вектор x, такой, что Ax = b
Обращение матрицы:
найти матрицу A⁻¹, такую, что AA⁻¹ = I

Для численных методов задача осложняется: даны значения $\hat{A}, \hat{b}$ с погрешностью, т.к. на практике не все действительные числа могут быть точно представлены в виде чисел с плавающей точкой; арифметические операции тоже вносят погрешность.

Векторные и матричные нормы

Векторная норма ‖·‖ : 𝐕 → ℝ_+:

Разделяет точки: ‖x‖ = 0 ⇔ x = 0
Однородна: λ ∈ ℂ ⇒ ‖λx‖ = |λ|·‖x‖
Удовлетворяет неравенству треугольника: ‖x + y‖ ⩽ ‖x‖ + ‖y‖

Матричная норма — п.1–3 и:

Удовлетворяет неравенству треугольника по произведению: ‖AB‖⩽‖A‖·‖B‖
(получается, что не всякая векторная норма является матричной)

Свойства:

‖I‖⩾1
‖A⁻¹‖·‖A‖⩾1

Подчинённые матричные нормы

Для нормы ‖ · ‖ на ℂⁿ определим матричную (подчинённую, индуцированную) норму на 𝐌_n по формуле:

‖A‖≡max_‖x‖=1‖Ax‖.

Для подчинённой нормы ‖I‖=1 и ‖Ax‖⩽‖A‖·‖x‖

Основные нормы для этого курса:

Максимальная столбцовая норма (𝐋¹): $‖A‖_{1} = \max\limits_{1⩽j⩽n} \sum\limits_{i=1}^n |a_{ij}| \qquad ‖x‖_{1} = Σ|x_{i}|$
Максимальная строчная норма (𝐋^∞): $‖A‖_{∞} = \max\limits_i \sum\limits_j |a_{ij}| \qquad ‖x‖_∞ = \max |x_{i}|$
Спектральная норма (𝐋²) (λ — собственное значение):

$$‖A‖_{2} = \max \{\sqrt{λ} \; | \; ∃v : Av = λv\} \qquad ‖x‖_2 = \sqrt{Σx_{i}^2}$$

Спектральная норма в первой задаче не используется.

Вычислительные ошибки и число обусловленности

Пусть x — точное решение системы Ax = b, а $\hat{x}$ — приближённое, полученное с помощью численных методов.
Назовём числом обусловленности (κ(A), cond(A)): $κ(A) ≡ ‖A‖·‖A^{-1}‖ \; \mbox{(для невырожденных матриц)} $

Невязка: $r ≡ b - A\hat{x}$

Относительная погрешность приближённого решения: $ \dfrac{‖x-\hat{x}‖}{‖x‖} ⩽ κ(A)\dfrac{‖r‖}{‖b‖}. $

Свойства числа обусловленности:

κ(A)⩾1
κ(A)=κ(A⁻¹)
κ(AB)⩽κ(A)κ(B)
$κ(A) ⩾ \left| \dfrac{λ_{max}(A)}{λ_{min}(A)} \right|$

Умножение матриц

Схема умножения

Умножение квадратных матриц

Тривиальная реализация:

for (i = 0; i < n; i++) {
  for (j = 0; j < n; j++) {
    c[i*n+j] = 0.;
    for (k = 0; k < n; k++) {
      c[i*n+j] = a[i*n+k] * b[k*n+j];

Блочная реализация (размер блока — порядка $\sqrt{\dfrac{L1}{3}}$)

for (i = 0; i < n/BLOCK; i++) {
  for (j = 0; j < n/BLOCK; j++) {
    zero_block(c, i, j);
    for (k = 0; k < n/BLOCK; k++) {
      multiply_block(a, b, temp, i, j, k);
      put_block(c, temp, i, j);

Умножение матриц (современное состояние исследований)

Вычислительная сложность тривиального алгоритма — O(n³) (т.е. если реализация умножения работает 1 с. на матрице размера 1000 × 1000, то она будет работать 8 с. на матрице размера 2000 × 2000, 64 с. на матрице размера 4000 × 4000 и т.д.). Развёртывание цикла, блочные оптимизации и т.п. могут только уменьшить константу.

Минимальный теоретически возможный порядок — O(n²), т.к. результат зависит от всех 2n² элементов входных матриц, но вопрос минимальной достижимой сложности на настоящее время не решён.

Лучший применимый на практике алгоритм — Strassen (1969), сложность — O(n^2.807355), есть проблемы с численной устойчивостью.

Лучший алгоритм — Coppersmith-Winograd, с последующими улучшениями (2014) сложность — O(n^2.3728639), но константа и производительность используемых операций не позволяет использовать его на практике с ощутимым преимуществом (проблемы с устойчивостью тоже никуда не деваются).

Невязка

Невязка (residual) — вектор (матрица) погрешности вычислений. Нас интересует норма невязки.

решение линейной системы: ‖Ax − b‖
умножить сгенерированную матрицу на какой-то*, тоже сгенерированный, вектор и вычесть правую часть (обычно будет один из векторов стандартного базиса — 0, …, 0, 1, 0, …, 0
обращение матрицы: ‖AA⁻¹ − I‖
умножить сгенерированную матрицу на другую* сгенерированную матрицу и вычесть единичную

Для «хороших» матриц наша цель по невязке порядка 10⁻¹⁶

Задание

Изменить программу, выполненную по заданию с прошлого занятия, следующим образом:

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

Список задач
или lynx в терминале.

Формулы для матриц

i, j = 0…(n − 1):

Симметричная: a_ij = |i − j| — по возможности тестируйте на больших размерах этой или следующей матрицы
Симметричная без нулей на диагонали: a_ij = |i − j|+1
Матрица Гильберта: $H_{ij} = \dfrac{1}{i+j+1}$ — нет смысла тестировать на размерах >100

Для нормы 𝐋² число обусловленности матрицы Гильберта 5 × 5:

κ(H)≈4.8 · 10⁵

Сложность алгоритма O(n^k) означает, что для задачи «размера» n выполнение алгоритма требует не более, чем C · n^k (члены младших степеней можно отбросить путём увеличения константы).