Ruby/Матрицы и векторы

Матрицы и векторы

править

Моё знакомство с реализацией матричного исчисления в Ruby началось ещё в мою бытность студентом МГИУ. И надо сказать, что это знакомство было неудачным. Дело в том, что простейшая программа вычисления определителя матрицы

require 'matrix'
p Matrix[[1, -2, 3], [3, 4, -5], [2, 4, 1]].det    #=> -50

( версия ruby1.9.1 выдает результат #=> (62/1) - прим. Sharipov.ru)

давала неверный результат (правильный ответ — 62). Как выяснилось позднее, эта проблема связана со спецификой целочисленной арифметики в Ruby (одна вторая в Ruby — нуль). Предположив это, я решил, что проблема легко решится, если я преобразую элементы матрицы к типу чисел с плавающей запятой (классу Float):

require 'matrix'
p Matrix[[1.0, -2.0, 3.0], [3.0, 4.0, -5.0], [2.0, 4.0, 1.0]].det    #=> 62.0

И в некоторых случаях это меня действительно выручало. Но не всегда. Стоило появиться делению на 3, как появлялись ненужные остатки и погрешности. И чем больше исходная матрица, тем больше вероятность появления таких остатков. В некоторых случаях это было несущественным, а в остальных приходилось прибегать к услугам специальных математических пакетов (например, Maxima). Было жутко неудобно и обидно за такую «кривую реализацию». (Я тогда писал курсовой, который решал все варианты контрольной работы по предмету «Математическое программирование». Да простит меня преподаватель, который так и не понял секрета тотальной успеваемости группы.)

На этом бы история и закончилась (как позже я узнал, на этом она заканчивалась для многих), но мне в руки попалась книга Programming Ruby 2ed с описанием возможностей стандартной библиотеки версии 1.8.2. Именно там (на странице 671) я наткнулся на описание библиотеки mathn. Уникальность её состоит в том, что она существенно расширяет возможности стандартных чисел, добавляя к ним рациональные числа (класс Rational) и комплексные числа (класс Complex).

Проще говоря, появляется возможность делить числа без погрешностей (класс Rational) и возможность извлекать квадратный корень из отрицательного числа (класс Complex).

Одновременно с этим она добавляет матричное и векторное исчисления (правда, почему-то в книге дополнительно подключали complex и matrix). И после этого матричное счисление начинает работать «из коробки», и ещё как работать! Хотите обратную матрицу? Пожалуйста! Хотите определитель? Нет ничего проще! Обидно только одно: программу к тому времени я уже написал и эти возможности мне были не нужны.

Чуть позднее, один из моих студентов написал мне письмо с просьбой объяснить как «работать с матрицами в Ruby». При этом он задал всего три вопроса:

  1. Как создать новую матрицу?
  2. Как послать в матрицу двумерный массив?
  3. Как поменять элемент матрицы?

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

Как создать новую матрицу?

править

Самый простейший способ создать матрицу — использовать метод «батарейка» (метод [] выглядит как индикатор заряда батареи на сотовом телефоне):

require 'mathn'
Matrix[[1, -2, 3], [3, 4, -5], [2, 4, 1]]

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

Как послать в матрицу двумерный массив?

править

Немножко перефразируем вопрос: как преобразовать двумерный массив в матрицу? Пусть ответ будет неожиданным, но это делается при помощи того-же метода «батарейка»:

require 'mathn'
array = [[1, -2, 3], [3, 4, -5], [2, 4, 1]]
Matrix[*array]

Оператор * преобразует array в список параметров, что позволяет существенно упросить процесс создания матрицы. Более того, это единственно удобный метод её создания. Все остальные методы создания матрицы (например, метод .new) работают не столь изящно и интуитивно.

Как изменить элемент матрицы?

править

И вот тут всплывает «недоделанность» матриц. Метода для изменения элемента матрицы в них нет! Для того, чтобы изменить элемент матрицы, надо преобразовать матрицу в массив, изменить элемент массива и преобразовать массив в матрицу. Примерно вот так (меняем элемент с индексом 0, 0):

require 'mathn'
matrix =  Matrix[[1, -2, 3], [3, 4, -5], [2, 4, 1]]
array = matrix.to_a
array[0][0] = 0
p Matrix[*array]    #=> Matrix[[0, -2, 3], [3, 4, -5], [2, 4, 1]]

Исправляем оплошность разработчиков

править

Для начала, рассматриваем поближе библиотеку matrix (исходник или описание в fxri) и выясняем, что для получения значения элемента используется метод «батарейка» с двумя аргументами. Вполне закономерно использовать метод «батарейка равно» (то есть []=) для изменения элемента. Сейчас мы его и реализуем:

require 'mathn'
class Matrix
    def []=(i, j, value)
        @rows[i][j] = value
    end
end

matrix =  Matrix[[1, -2, 3], [3, 4, -5], [2, 4, 1]]
matrix[0, 0] = 5
p matrix  #=> Matrix[[5, -2, 3], [3, 4, -5], [2, 4, 1]]

Почему не могли этого сделать разработчики, я так и не понял. Скорее всего по идеологическим соображениям («не дело, чтобы матрицы вели себя как простые массивы»).
в версиях 1.9.x и выше последние 3 строчки работают, без определения класса. posix_nik
в версии 1.9.2-p136 метод []= почему-то приватный

Методы работы с матрицами

править

Из методов работы с матрицами наиболее важные — вычисление определителя (.det), вычисление обратной матрицы (.inv), поиск минора матрицы (.minor), преобразование матрицы в массив (.to_a), умножение (*), сложение (+), вычитание (-) и деление (/) матриц. Без остальных методов можно обойтись, поэтому остальные методы изучите самостоятельно.

Зачем нужны векторы?

править

Во-первых, для Ruby вектор — это объект класса Vector. Подключается он одновременно с матрицами (класс Matrix). Во-вторых, вектор очень похож на массив, но с одним существенным отличием (определяющем полезность вектора): массивы и векторы по-разному складываются и вычитаются. Давайте рассмотрим небольшой пример:

require 'mathn'
array = [1, 2, 3, 4, 5]
vector = Vector[*array]
p vector + vector    #=> Vector[2, 4, 6, 8, 10]
p array + array      #=> [1, 2, 3, 4, 5, 1, 2, 3, 4, 5]

Обратите внимание, что Vector создаётся точно также, как и матрица. По сути, вектор — это матрица, которая состоит лишь из одной строки. А матрица, в свою очередь, — это массив векторов.

В чём разница? Для сложения векторов необходимо сложить их соответствующие координаты. В случае же с массивами происходит конкатенация (сцепление) массивов.

Хитрости

править

Задачи

править

Решение систем линейных уравнений методом Гаусса

править

Реализуем метод Гаусса.

require 'mathn'
equation = [Vector[1, 2, 1, 1], Vector[1, 5, 6, 2], Vector[1, 5, 7, 10]]

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

require 'mathn'
equation = [[1, 2, 1, 1], [1, 5, 6, 2], [1, 5, 7, 10]].map{ |array| Vector[*array] }

# Прямой ход метода Гаусса
equation[0] /= equation[0][0]

И вот тут нас ждет неприятный сюрприз: у векторов не реализован метод /. Ну и ладно, мы, программисты, не гордые. Сами напишем!

require 'mathn'

class Vector
    def /(arg)
        self * (1 / arg)
    end
end

equation = [[1, 2, 1, 1], [1, 5, 6, 2], [1, 5, 7, 10]].map{ |array| Vector[*array] }

# Прямой ход метода Гаусса
equation[0] /= equation[0][0]

Теперь работает (только со скалярами, но нам больше и не надо). Заканчиваем реализовывать прямой проход метода Гаусса:

require 'mathn'

class Vector
  def /(arg)
    self * (1/arg)
  end
end

equation = [[1, 2, 1, 1], [1, 5, 6, 2], [1, 5, 7, 10]].map{ |array| Vector[*array] }

# Прямой ход метода Гаусса
equation[0] /= equation[0][0]
equation[1] -= equation[0] * equation[1][0]
equation[2] -= equation[0] * equation[2][0]

equation[1] /= equation[1][1]
equation[2] -= equation[1] * equation[2][1]

equation[2] /= equation[2][2]

p equation  #=> [Vector[1, 2, 1, 1], Vector[0, 1, 5/3, 1/3], Vector[0, 0, 1, 8]]

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

require 'mathn'

class Vector
  def /(arg)
    self * (1/arg)
  end
end

equation = [[1, 2, 1, 1], [1, 5, 6, 2], [1, 5, 7, 10]].map{ |array| Vector[*array] }

# Прямой ход метода Гаусса
equation[0] /= equation[0][0]
equation[1] -= equation[0] * equation[1][0]
equation[2] -= equation[0] * equation[2][0]

equation[1] /= equation[1][1]
equation[2] -= equation[1] * equation[2][1]

equation[2] /= equation[2][2]

# Обратный ход метода Гаусса
equation[1] -= equation[2] * equation[1][2]
equation[0] -= equation[2] * equation[0][2]

equation[0] -= equation[1] * equation[0][1]

p equation  #=> [Vector[1, 0, 0, 19], Vector[0, 1, 0, -13], Vector[0, 0, 1, 8]]

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

p equation.map{ |vector| vector.to_a }.transpose[-1]  #=> [19, -13, 8]

Теперь задачу можно считать решённой. Жаль только, что программа работает только для уравнения 3×4. Надо бы добавить несколько итераторов, чтобы они самостоятельно определяли размерность уравнения. Для этого нужно проследить чередование индексов и записать для них итераторы:

require 'mathn'

class Vector
    def /(arg)
        self * (1/arg)
    end
end

equation = Matrix[[1, 2, 1, 1], [1, 5, 6, 2], [1, 5, 7, 10]].row_vectors

# Прямой ход метода Гаусса
(0...equation.size).each{ |i|
    equation[i] /= equation[i][i]
    (i+1...equation.size).each{ |j| equation[j] -= equation[i] * equation[j][i] }
}

# Обратный ход метода Гаусса
(1...equation.size).to_a.reverse.each{ |i|
    (0...i).each{ |j| equation[j] -= equation[i] * equation[j][i] }
}

p equation.map{ |vector| vector[-1] }  #=> [19, -13, 8]

Обратите внимание, что equation задаётся через матрицу (которая, не без помощи метода .row_vectors, преобразуется в массив векторов). Также обратите внимание, что получить последний столбец можно посредством итератора .map и метода «батарейка».