Библиотека сайта rus-linux.net
Глава 11
Numeric Python
Numeric Python (численный Python) - собрание модулей, в центре которого быстродействующее
ядро для работы с многомерными массивами, необходимое для многих численных приложений.
В какой-то степени Numeric вносит в Python возможности таких пакетов и систем, как
MatLab, Octave (аналог MatLab), APL, J, S+, IDL. Пользователи этих пакетов без труда
обнаружат в Numeric знакомые черты.
Numeric Python имеет средства для:
?матричных вычислений (LinearAlgebra)
?быстрого преобразования Фурье (FFT)
?работы с недостающими экспериментальными данными (MA)
?статистического моделирования (RNG)
В необходимых случаях Numeric может быть сконфигурирован для работы со ставшими
классическими численными библиотеками LAPACK и FFTPACK
(в их отсутствие используется облегченная версия LAPACK - LALITE). Кроме того,
имеется специальный модуль MLab для эмуляции базовых функций программы MatLab.
Ниже будут коротко рассмотрены возможности модулей Numeric и LinearAlgebra.
11.1. Модуль Numeric
Модуль Numeric вносит в Python полноценный объект-массив и большое число функций
для оперирования с ним.
Массив - набор однородных элементов. По сравнению с массивами модуля array, массивы модуля
Numeric могут иметь более одной размерности.
Массив можно создать с помощью функции array():
from Numeric import *
print array([[1, 2], [3, 4], [5, 6]])
print array([[1, 2, 3], [4, 5, 6]], Float)
print array([65, 66, 67, 68, 69], 'c')
Результат выполнения:
[[1 2]
[3 4]
[5 6]]
[[ 1. 2. 3.]
[ 4. 5. 6.]]
[A B C D E]
Для элементов массивов можно использовать типы Int8-Int32, UnsignedInt8, Float8-Float64,
Complex8-Complex64 и PyObject. Числа 8, 16, 32 и 64 обозначают количество бит для
хранения величины. Типы Int, Float и Complex соответствуют наибольшим принятым на
данной платформе значениям. В массивах можно также хранить ссылки на произвольные
объекты (тип PyObject), но это, естественно, будет сказываться на производительности.
Для некоторых применений операции, существующие для числовых массивов, очень удобны и
для массивов с произвольными объектами.
Количество размерностей (осей) и длина массива по каждой оси называются формой (shape)
массива. Прочитать или задать форму можно с помощью атрибута shape:
from Numeric import *
a = array(range(12), Int)
print "a.shape:", a.shape
print a
a.shape = (3, 4)
print "a.shape:", a.shape
print a
Результат выполнения:
a.shape: (12,)
[ 0 1 2 3 4 5 6 7 8 9 10 11]
a.shape: (3, 4)
[[ 0 1 2 3]
[ 4 5 6 7]
[ 8 9 10 11]]
Для задания формы можно использовать и функцию Numeric.reshape(), причем, действовать
она может сразу над последовательностями, создавая объект-массив.
>>> import Numeric
>>> print Numeric.reshape("правда", (2, -1))
[[п р а]
[в д а]]
В этом примере -1 в указании формы - не случайность. Общее количество элементов массива
известно (в примере - 6), поэтому длину одной (и не более одной!) из осей можно
не задавать, записывая на ее месте -1.
Атрибут flat выпрямляет массив в одномерный:
>>> import Numeric
>>> a = array([[1, 2], [3, 4]])
>>> print a.flat
[1 2 3 4]
Похожа на Numeric.reshape() функция Numeric.resize() для придания массиву необходимой
формы и размера:
>>> import Numeric
>>> print Numeric.resize("правда", (10,))
[п р а в д а п р а в]
Как уже говорилось при рассмотрении объектов-срезов, объекты-массивы Numeric используют
расширенный синтаксис выделения среза. Следующие наборы примеров хорошо иллюстрируют
различные виды записи срезов. Здесь же мы видим функцию Numeric.arrayrange(), которая
является аналогом range() для массивов.
import Numeric
a = Numeric.arrayrange(24) + 1 # массив чисел от 1 до 24
a.shape = (4, 6) # форма 4x6
examples = [ ("a[1,2]", "элемент 1,2"),
("a[1,:]", "строка 1"),
("a[1]", "тоже строка 1"),
("a[:,1]", "столбец 1"),
("a[-2,:]", "предпоследняя строка"),
("a[0:2,1:3]", "окно 2x2"),
("a[1,::3]", "каждый третий элемент из строки 1"),
("a[:,::-1]", "реверс элементов в каждой строке"),
]
print a
for (expr, expl) in examples:
print "%s (%s):\n%s" % (expr, expl, str(eval(expr)))
Результат выполнения:
[[ 1 2 3 4 5 6]
[ 7 8 9 10 11 12]
[13 14 15 16 17 18]
[19 20 21 22 23 24]]
a[1,2] (элемент 1,2):
9
a[1,:] (строка 1):
[ 7 8 9 10 11 12]
a[1] (тоже строка 1):
[ 7 8 9 10 11 12]
a[:,1] (столбец 1):
[ 2 8 14 20]
a[-2,:] (предпоследняя строка):
[13 14 15 16 17 18]
a[0:2,1:3] (окно 2x2):
[[2 3]
[8 9]]
a[1,::3] (каждый третий элемент из строки 1):
[ 7 10]
a[:,::-1] (реверс элементов в каждой строке):
[[ 6 5 4 3 2 1]
[12 11 10 9 8 7]
[18 17 16 15 14 13]
[24 23 22 21 20 19]]
Второй набор примеров использует многоточие (Ellipsis). Многоточие ставится для указания
произвольного числа осей, т. е. вместо :,:,...,:.
import Numeric
a = Numeric.arrayrange(24) + 1 # массив чисел от 1 до 24
a.shape = (2,2,2,3) # форма 2x2x2x3
examples = [ ("a[0,...]", "0-й блок"),
("a[0,:,:,0]", "срез по первой и последней осям"),
("a[0,...,0]", "то же"),
]
print a
for (expr, expl) in examples:
print "%s (%s):\n%s" % (expr, expl, str(eval(expr)))
Результат выполнения:
[[[[ 1 2 3]
[ 4 5 6]]
[[ 7 8 9]
[10 11 12]]]
[[[13 14 15]
[16 17 18]]
[[19 20 21]
[22 23 24]]]]
a[0,...] (0-й блок):
[[[ 1 2 3]
[ 4 5 6]]
[[ 7 8 9]
[10 11 12]]]
a[0,:,:,0] (срез по первой и последней осям):
[[ 1 4]
[ 7 10]]
a[0,...,0] (то же):
[[ 1 4]
[ 7 10]]
Для порождения массивов, состоящих из одних нулей или единиц, можно использовать функции
zeros() и ones() соответственно. Для порождения единичной матрицы размера n можно
воспользоваться identity(n).
Отдельно требуется затронуть вопрос о копировании массивов. Дело в том, что во
многих операциях, которые только меняют форму или выделяют некоторые срезы массива,
результат представляет лишь еще одно "окно" к тем же данным, и изменение данных в
этом окне сказывается на исходном массиве. Следующий пример показывает, что для
копирования массивов лучше явно использовать copy() из модуля copy.
from Numeric import *
from copy import copy
a = array([1, 2, 3, 4, 5, 6, 7, 8])
print a
b, c, d = a[:2], reshape(a, (4, 2)), copy(a)
b[0], c[0, 1], d[2] = -1, -2, -3
print a
Результат выполнения:
[1 2 3 4 5 6 7 8]
[-1 -2 3 4 5 6 7 8]
Функция array() для создания массива также делает копию, если ее аргумент -
массив. Кроме того, имеется функция asarray(), которая не создает нового массива,
когда ее аргумент уже массив.
Модуль Numeric определяет целый набор функций, аналогичных соответствующим встроенным
функциям, функциям из модуля math и др. В табл. 11.1 приведен список доступных функций,
называемых универсальными в силу того, что они применимы не только к массивам, но и к
последовательностям (в результате все равно получается массив). См. также разд. 9.4.1,
где перечислены функции модуля math.
Таблица 11.1. Универсальные функции модуля Numeric
Название функции
Действие функции
add(), subtract()
Сложение, вычитание
multiply(), divide(), remainder(), fmod()
Умножение, деление, получение остатка от деления (для целых чисел и чисел с плавающей
запятой)
power()
Возведение в степень
sqrt()
Корень квадратный
negative(), absolute(), fabs()
Смена знака и абсолютное значение
ceil(), floor()
Наименьшее (наибольшее) целое, большее (меньшее) или равное аргументу
hypot()
Длина гипотенузы (даны два катета)
sin(), cos(), tan()
Тригонометрические функции
arcsin(), arccos(), arctan()
Обратные тригонометрические функции
arctan2()
Арктангенс от частного аргументов
sinh(), cosh(), tanh()
Гиперболические функции
arcsinh(), arccosh(), arctanh()
Обратные гиперболические функции
exp()
Экспонента (e в степени)
log(), log10()
Натуральный и десятичный логарифмы
maximum(), minimum()
Максимум и минимум
conjugate()
Сопряжение (для комплексных чисел)
equal(), not_equal()
Равно, не равно
greater(), greater_equal()
Больше, больше или равно
less(), less_equal()
Меньше, меньше или равно
logical_and(), logical_or(), logical_xor(), logical_not()
Логические И, ИЛИ, исключающее ИЛИ и отрицание
Таблица 11.1 (окончание)
Название функции
Действие функции
bitwise_and(), bitwise_or(), bitwise_xor(), invert()
Побитовые И, ИЛИ, исключающее ИЛИ и инверсия
left_shift(), right_shift()
Побитовые сдвиги влево и вправо
Эти функции являются объектами типа ufunc и применяются к массивам поэлементно. Кроме
того, их специальные методы accumulate, outer, reduce и reduceat предназначены для
выполнения операций над массивами в определенном режиме (соответственно): аккумулирование
результата, внешнее "произведение", сокращение и сокращение с заданными точками.
В представленных функциях и их применении над массивами в различных режимах кроется
огромный потенциал. С помощью умелого подбора нужной комбинации операций можно
запрограммировать весьма нетривиальные алгоритмы, которые будут выполняться гораздо
быстрее, чем их аналоги, записанные с применением вложенных циклов. Однако проще
всего объяснить действие функции на примере сложения add.
>>> from Numeric import add
>>> add([[1, 2],
... [3, 4]], [[1, 0],
... [0, 1]])
[[2 2]
[3 5]]
>>> add([[1, 2],
... [3, 4]], [1, 0])
[[2 2]
[4 4]]
>>> add([[1, 2], [3, 4]], 1)
[[2 3]
[4 5]]
>>> add.reduce([1,2,3,4]) # т. е. 1+2+3+4
10
>>> add.reduce([[1, 2], # т. е. [1+3 2+4]
... [3, 4]], 0)
[4 6]
>>> add.reduce([[1, 2], # т. е. [1+2 3+4]
... [3, 4]], 1)
[3 7]
>>> add.accumulate([1,2,3,4]) # т. е. [1 1+2 1+2+3 1+2+3+4]
[ 1 3 6 10]
>>> add.reduceat(range(10),[0,3,6]) # т. е. [0+1+2 3+4+5 6+7+8+9]
[ 3 12 30]
>>> add.outer([1,2], [3,4]) # т. е. [[1+3 1+4] [2+3 2+4]]
[[4 5]
[5 6]]
Методы accumulate(), reduce() и reduceat() принимают третий, необязательный, аргумент
- номер оси, которая используется для соответствующего действия. По умолчанию
используется нулевая ось.
Сами универсальные функции помимо одного или двух необходимых аргументов позволяют
задавать и еще один аргумент, который принимает результат производимой функцией
операции. Тип третьего аргумента должен строго соответствовать типу результата.
Например sin, взятый от целых чисел, может иметь тип Float.
>>> from Numeric import array, sin, Float
>>> a = array([0, 1, 2])
>>> r = array([0, 0, 0], Float)
>>> sin(a, r) # sin(a) поместить в r
>>> print r
[ 0. 0.84147098 0.90929743]
В табл. 11.2 приводятся функции модуля Numeric, которые являются краткой записью
некоторых наиболее употребительных сочетаний функций и методов из рассмотренных
выше.
Таблица 11.2. Некоторые функции модуля Numeric и их аналоги
Функция и ее аргументы
Аналог функции
sum(a, axis)
add.reduce(a, axis)
cumsum(a, axis)
add.accumulate(a, axis)
product(a, axis)
multiply.reduce(a, axis)
cumproduct(a, axis)
multiply.accumulate(a, axis)
alltrue(a, axis)
logical_and.reduce(a, axis)
sometrue(a, axis)
logical_or.reduce(a, axis)
Следующий пример показывает, что всегда следует внимательно смотреть, где размещается
результат операции:
Листинг 11.1. Клеточный автомат
from Numeric import *
lenta, lenta[5] = zeros(16), 1
print "Пытаемся делать сложение на месте:"
for i in xrange(3):
print lenta
add(lenta[:-1], lenta[1:], lenta[1:])
lenta, lenta[5] = zeros(16), 1
print "Теперь используем промежуточный массив:"
for i in xrange(3):
print lenta
lenta[1:] = add(lenta[:-1], lenta[1:])
lenta, lenta[5] = zeros(16), 1
a = zeros(15)
print "Используем явный промежуточный массив:"
for i in xrange(3):
print lenta
add(lenta[:-1], lenta[1:], a)
lenta[1:] = a
Результат выполнения:
Пытаемся делать сложение на месте:
[0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1]
[ 0 0 0 0 0 1 2 3 4 5 6 7 8 9 10 11]
Теперь используем промежуточный массив:
[0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0]
Используем явный промежуточный массив:
[0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0]
В первом случае мы попытались сэкономить на памяти и получили испорченный результат.
Второй и третий варианты избежали ошибки, т. к. результат операции в них сохраняется
в промежуточном массиве. К счастью, в нашем случае ошибку легко увидеть, чего не
скажешь о сложных матричных алгоритмах.
Для сохранения массивов (и других объектов) в файлоподобных объектах можно либо
использовать функции Numeric.dump() и Numeric.load(), либо универсальные функции
модуля pickle (см. разд. 9.2.7). Судя по исходным текстам модуля Numeric, лучше
воспользоваться первым вариантом:
from Numeric import *
file = open('my.dat', 'wb')
a1 = array(range(500), Float)
a2 = reshape(range(250), (25, -1))
a3 = {'my' : 'data'}
dump((a1, a2, a3), file)
print a1[100], a2[0, 0], a3
file.close()
# возможно, уже в другой программе, даже на другой машине:
file = open('my.dat', 'rb')
(a, b, c) = load(file)
print a[100], b[0, 0], c
Настал черед рассказать о других функциях работы с массивами.
Функция Numeric.take() позволяет (дословно) взять часть массива по заданным на
определенной оси индексам. По умолчанию номер оси (третий аргумент) равен 0.
>>> from Numeric import *
>>> a = reshape(arrayrange(25), (5, 5))
>>> print a
[[ 0 1 2 3 4]
[ 5 6 7 8 9]
[10 11 12 13 14]
[15 16 17 18 19]
[20 21 22 23 24]]
>>> print take(a, [1], 0)
[ [5 6 7 8 9]]
>>> print take(a, [1], 1)
[[ 1]
[ 6]
[11]
[16]
[21]]
>>> print take(a, [[1,2],[3,4]])
[[[ 5 6 7 8 9]
[10 11 12 13 14]]
[[15 16 17 18 19]
[20 21 22 23 24]]]
В отличие от операции выделения среза, функция Numeric.take() оставляет размерность
массива той же, если конечно, структура заданных индексов одномерна. Результат
take(a, [[1,2],[3,4]]) показывает, что взятые по индексам части помещаются в
массив со структурой самих индексов, как если вместо 1 написать [ 5 6 7 8 9],
вместо 2 взять [10 11 12 13 14] и т. д.
Функция Numeric.diagonal() возвращает диагональ.
diagonal(a[, смещ[, ось1[, ось2]]]) -> a1
где a - исходный массив; смещ - смещение вправо от "главной" диагонали (по умолчанию
0); ось1 - первая из осей, на которых берется диагональ (по умолчанию 0); ось2 -
вторая ось (по умолчанию 1); a1 - массив с результатом.
>>> from Numeric import *
>>> a = reshape(arrayrange(16), (4, 4))
>>> print a
[[ 0 1 2 3]
[ 4 5 6 7]
[ 8 9 10 11]
[12 13 14 15]]
>>> for i in range(-3, 4):
... print diagonal(a, i)
...
[12]
[ 8 13]
[ 4 9 14]
[ 0 5 10 15]
[ 1 6 11]
[2 7]
[3]
Наверное, к этому моменту идея модуля Numeric стала понятна. В табл. 11.3 приведены
описания функций модуля.
Таблица 11.3. Функции модуля Numeric
Функция и ее аргументы
Назначение функции
array(a[, тип
[, savespace]])
Создание массива на основе последовательности a данного типа. Аргумент savespace
является флагом, указывающим необходимость сохранения места
zeros(форма[, тип])
Массив нулей заданной формы и типа
ones(форма[, тип])
Массив единиц заданной формы и типа
fromstring(строка
[, колво[, тип]])
Создание массива на основе бинарных данных, хранящихся в строке
take(a, индексы[, ось])
Выбор частей массива a на основе индексов по оси
put(a, индексы, b)
Присваивание частям массива, a[n] = b[n] для всех индексов n
putmask(a, маска, b)
Присваивание a элементов из b, для которых маска имеет значение "истина"
reshape(a, форма)
Возвращает массив нужной формы (нового массива не создает)
repeat(a, n[, ось])
Повторяет элементы массива a n раз по оси
choose(a, (b0,...,bn))
Конструирует массив на основе элементов, взятых по индексам из a (индексы от 0 до
n включительно). Формы массивов a, b1, ..., bn должны совпадать
cross_correlate(a, b
[, режим])
Взаимная корреляция двух массивов. Параметр режим может быть равен 0, 1 или 2
searchsorted(a, i)
Для каждого элемента из i найти его место в массиве a. Массив a должен быть
одномерным и отсортированным. Результат имеет форму массива i
sum(a[, ось])
Суммирование по оси ось массива a
cumsum(a[, ось])
Суммирование с промежуточными результатами
product(a[, ось])
Произведение по оси ось массива a
cumproduct(a[, ось])
Произведение по оси ось массива a с промежуточными результатами
alltrue(a[, ось])
Логическое И по всей оси ось массива a
sometrue(a[, ось])
Логическое ИЛИ по всей оси ось массива a
allclose(a, b
[, eps[, A]])
Сравнение a и b с заданными относительными eps и абсолютными A погрешностями.
По умолчанию eps равен 1.0e-1, а A - 1.0e-8
Таблица 11.3 (продолжение)
Функция и ее аргументы
Назначение функции
arrayrange(start
[, stop[, step[, тип]]]) или arange()
Аналог range() для массивов
asarray(a[, тип
[, savespace]])
То же, что и array, но не создает новый массив, если a - уже массив
convolve(a, b[, режим])
Свертка двух массивов. Аргумент режим может быть равен 0, 1 или 2
swapaxes(a, ось1, ось1)
Смена осей (частный случай транспонирования)
concatenate(a[, ось])
Соединение двух массивов (конкатенация) по заданной оси (по умолчанию - по
нулевой оси).
transpose(a[, оси])
Перестановка осей в соответствии с аргументом оси, либо, если оси не заданы -
расположение их в обратном порядке
sort(a[, ось])
Сортировка элементов массива по заданной оси
argsort(a[, ось])
Индексы отсортированного массива, такие, что take(a,argsort(a, axis),axis) дает
отсортированный массив a, как если бы делалось sort(a, axis)
argmax(a[, ось])
Индекс максимального значения в массиве на заданной оси
argmin(a[, ось])
Индекс минимального значения в массиве на заданной оси
innerproduct(a, b)
Внутреннее произведение двух массивов (по общему измерению). Для успеха операции
a.shape[-1] должен быть равен b.shape[-1]. Форма результата будет равна
a.shape[:-1] + b.shape[:-1]. Элементы пропадающего измерения попарно умножаются
и суммируются
dot(a, b)
Внутреннее (матричное) произведение массивов. По определению: innerproduct(a,
swapaxes(b, -1, -2)), т. е. с переставленными последними осями, как и подобает
произведению матриц
outerproduct(a, b)
Внешнее произведение a и b
resize(a, форма)
Возвращает массив с произвольной новой формой
indices(измерения
[, тип])
Возвращает массив сеток индексов заданной длины по каждому измерению с изменением
поочередно по каждому измерению. Например, indices([2, 2])[1] дает двумерный
массив [[0, 1], [0, 1]]
Таблица 11.3 (продолжение)
Функция и ее аргументы
Назначение функции
fromfunction(f,
измерения)
Строит массив путем вызова функции f(), в качестве аргументов которой выступают
значения кортежа сеток индексов. Фактически является сокращением
для apply(f, tuple(indices
(измерения)))
diagonal(a[, k[, ось1
[, ось2]]])
Взятие k-й диагонали массива a в плоскости осей ось1 и ось2
trace(a[, k[, ось1[, ось2]]])
Сумма элементов вдоль диагонали. То есть add.reduce(diagonal(a, k, ось1, ось2))
dump(obj, file)
Запись массива a (в двоичном виде) в открытый файлоподобный объект file. Файл
должен быть открыт в бинарном режиме. Можно записать несколько объектов подряд
dumps(obj)
Строка с двоичным представлением объекта obj
load(file)
Чтение массива из файла (или файлоподобного объекта) file. Файл должен быть открыт
в бинарном режиме
loads(s)
Возвращает объект, соответствующий бинарному представлению в строке
ravel(a)
Распрямление массива в одномерный. Аналогично reshape(a, (-1,))
nonzero(a)
Возвращает индексы ненулевых элементов одномерного массива
shape(a)
Возвращает форму массива a
where(массив_условие, м1, м2)
Выбор элементов на основании массив_условие из м1 (если не ноль) и м2 (при нуле)
поэлементно. Равносилен choose(not_equal(condition, 0), (y, x)). Формы
массивов-аргументов должны совпадать
compress(условие,
a[, ось])
Возвращает массив только из тех элементов массива a, для которых условие
истинно (не ноль)
clip(a, a_min, a_max)
Обрубает значения массива a так, чтобы они находились между значениями из a_min
и a_max поэлементно. Если элемент выходит за пределы, указанные элементами с теми
же индексами в a_min и a_max, он принимает ближайшее значение из этого промежутка
zeros(форма[, тип])
Массив из нулей заданной формы и типа
Таблица 11.3 (окончание)
Функция и ее аргументы
Назначение функции
ones(форма[, тип])
Массив из единиц заданной формы и типа
identity(n)
Возвращает двумерный массив формы (n, n)
Модуль Numeric также определяет константы e (число e) и pi (число пи).
Конечно, научиться применять эти функции можно только в процессе практики, при
решении конкретной задачи.
11.2. Модуль LinearAlgebra
Модуль LinearAlgebra предоставляет доступ к некоторым часто использующимся в
линейной алгебре алгоритмам. Предполагается, что читатель хоть немного знаком
с линейной и матричной алгеброй и без труда найдет применение описываемым ниже
функциям.
Функция determinant() находит детерминант (определитель) матрицы:
>>> import Numeric, LinearAlgebra
>>> print LinearAlgebra.determinant(
... Numeric.array([[1, -1],
... [2, 3]]))
5
Функция solve_linear_equations() позволяет решать линейные уравнения вида ax=b
по данным ей аргументам a и b. Например:
>>> import Numeric, LinearAlgebra
>>> a = Numeric.array([[1.0, 2.0],
... [0.0, 1.0]])
>>> b = Numeric.array([1.1, 0.5])
>>> x = LinearAlgebra.solve_linear_equations(a, b)
>>> print "x =", x
x = [ 0.1 0.5]
>>> print "Проверка:", Numeric.dot(a, x) - b
Проверка: [ 0. 0.]
В случае, когда матрица a имеет нулевой определитель, возбуждается исключение
LinearAlgebraError с уточнением: Singular matrix.
Функция inverse() позволяет обращать матрицу. Только не следует решать линейные
уравнения с помощью LinearAlgebra умножением на обратную матрицу, ведь она сама
определена так:
def inverse(a):
return solve_linear_equations(a, Numeric.identity(a.shape[0]))
Функция eigenvalues() находит собственные значения матрицы, а eigenvectors() -
пару: собственные значения, собственные векторы. Следующий пример демонстрирует
эти две функции:
from Numeric import *
from LinearAlgebra import *
a = array([[-5, 2],
[2, -2]])
lmd = eigenvalues(a)
print "Собственные значения:", lmd
(lmd, v) = eigenvectors(a)
print "Собственные вектора:"
print v
print "Проверка:", dot(a, v[0]) - v[0] * lmd[0]
Результат выполнения:
Собственные значения: [-1. -6.]
Собственные вектора:
[[ 0.4472136 0.89442719]
[-0.89442719 0.4472136 ]]
Проверка: [ 6.66133815e-16 -3.33066907e-16]
В последней строке числа совсем маленькие (практически нули), значит, собственные
числа и вектора найдены верно.
Еще одна функция, использующаяся для исследования матриц - сингулярное разложение
матрицы.
singular_value_decomposition(A[, full]) -> (U, d, Vt)
где A - исходная матрица, которую представляют в виде U*diag(d)*Vt; full - влияет
на размерности результирующих матриц, если A не является квадратной, по умолчанию 0;
U - левый множитель (ортогональная матрица); d - вектор диагонали; Vt - правый
множитель (ортогональная матрица).
Рассмотрим пример:
Листинг 11.2. Сингулярное разложение
from Numeric import *
from LinearAlgebra import *
def diag(d):
"""Породить диагональную матрицу по заданной диагонали"""
n = d.shape[0]
r = zeros([n]*2, d.typecode())
ravel(r)[0::n+1] = d
return r
A = array([[1, 2],
[-1, 4]])
(U, d, Vt) = singular_value_decomposition(A)
print "Проверяем ортогональность матриц U и Vt:"
print abs(determinant(U)), abs(determinant(Vt))
print "Проверяем правильность сингулярного разложения:"
print dot(dot(U, diag(d)), Vt) - A
Результат выполнения:
Проверяем ортогональность матриц U и Vt:
1.0 1.0
Проверяем правильность сингулярного разложения:
[[ 0. 0.]
[ 0. 0.]]
Модуль LinearAlgebra является облегченной версией (эмуляцией) модуля LinAlg.
Последний же является высокоуровневым интерфейсом к хорошо отлаженной и
эффективной библиотеке LAPACK, написанной на FORTRAN или C.