Matlab

Материал из MachineLearning.

Matlab — язык программирования и система научных и инженерных расчетов, построенная на основе интерпретатора этого языка. Matlab, сокращение от «Matrix Laboratory», предназначен в первую очередь для выполнения алгоритмов, использующих векторы и матрицы.

Язык программирования Matlab (иногда также называется M-code) изначально был разработан с целью упрощения работы с процедурами широко распространенной в 70-80 годы библиотеки алгоритмов линейной алгебры LINPACK. Впоследствии он развился в мощный язык с богатым набором типов данных.

Оболочка Matlab состоит из командной строки, текстового редактора со встроенным отладчиком и окнами со списком файлов, списком видимых переменных и с историей введенных команд.

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

Краткий обзор языка Matlab

Векторы и матрицы

ввести вектор-столбец

vec  = [1; 2; 3]

вектор-строку

vec = [1 2 3]

и матрицу

mat  = [1 2 3; ...
        4 5 6]

задать матрицы

m = 3;
n = 2;
mat1 = ones(m,n)
mat0 = zeros(m,n)
mat2 = rand(n,n)
mat4 = randn(n,n)
mat3 = rand(n,n,n)

задать векторы

vec = [0:0.01:1]
vec = [0:100]
vec = linspace(0,1,100)

соединить векторы

vec = [0:0.5:1, 2:0.4:3]

добавть элемент к вектору-столбцу

vec(end)
vec(end+1) = 3

добавить элемент в вектору-строке

vec = [-1 vec]

векторизовать матрицу

vec = mat4(:)

удалить элементы

idx = find(vec < 0.05)
vec(idx) = []
vec=[]

соединить матрицы

mat1 = [1:5]'*[1:5]
mat2 = mat1(2:4,2:5)
mat2 = mat2'
mat2 = [mat1; mat1]
mat2 = [mat1, mat1]

текстовая строка - это вектор

i = 2
str = ['The lucky number is ', num2str(i)]

но не наоборот

str = '';              % empty string
str = sprintf(' 1 plus%s %d',str,i)

Операции

скалярное произведение

x = [1:5]';
norm = x'*x
table = x*x'

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

A = rand(5,2)
ATA = A'*A
ATA1 = inv(ATA)
ATA1'*ATA

поэлементные операции

y = 1+x
y = x+x
y = x.*x
Y = A.*A

Типы данных

структура

PP.niter = 100;
PP.tolX = 10e-6;
PP.display = 1;
niter = PP(1).niter;
PP

массив ячеек

strarray = {'one','two','et cetera'}
strarray = {'north','south';'west','east'}
fprintf(1, 'Go to the %s!\n', strarray{1,1});

возвращает ячейку

strarray(1,1)          % note that the following item is a cell

возвращает содержимое

strarray{1,1}

Цикл

цикл for

x=rand(5,1);				% get some x
for i=1:length(x)
	fprintf(1,'%1.2f\n',x(i));	% some code uses x(i)
end

итератор может быть элементом вектора

x=rand(5,1);				% get some x
idx = find(x>0.5);			% find some elements of x to process
for i = idx				% here idx is integer-valued vector
	fprintf(1,'%1.2f\n',x(i));	% some code uses x(i)
end

или вектором-столбцом матрицы

mat = [1 3 5 7; 2 4 6 8]
for i = mat				% no need to use additional lines if you deal with vectors
	fprintf(1,'i'' = %s\n',num2str(i'));
end

см. также

doc while % etc.

Условный переход

оператор if использует как логические так и целочисленные переменные; можно использовать операции 'and(a, b)', синоним 'a && b' и 'or(a, b)', синоноим 'a || b'

a = 1;
b = 2;
c = [];		% empty
d = NaN;		% Non a Number
e = -Inf;		% Infinity
if or (a>b, or( isempty(c), or( isnan(d), isinf(e))))
	disp('something remarkable happens');
end

но удобнее использовать

if any([a>b, isempty(c), isnan(d), isinf(e)])
	disp('something remarkable happens')
end

Эффективное программирование

Matlab в своих базовых конструкциях сложнее C/C++, поэтому желательно смотреть чужой код с целью выявления удачных приемов.

Пример 1. Примечание: представленные в обзоре времена работы зависят от конкретного компьютера, системы, на которой он работает, загруженности компьютера в момент вычисления примера, числа свободных положительных зарядов на шерсти Вашего кота... Короче, много от чего. Поэтому если при запуске скрипта на Вашем компьютере времена работы не будут совпадать с тем, что вы видете в этом обзоре - не пугайтесь, это нормально и ожидаемо. Главная мысль состоит в том, что правильно написанный код на MATLAB начинает не только быстрее работать, но и (почти всегда) лучше выглядеть и читаться.

Игрушечный пример для иллюстрации способов работы с функциями и векторами в среде MATLAB

size_of_vector = 10000;
random_vector = rand(1, size_of_vector);

Получили вектор, компоненты которого - случайные числа, полученные из равномерного на отрезке [0,1] распределения. Теперь возьмем от каждого его элемента логарифм различными способами.

Способ первый - совсем плохой, самый медленный - в цикле наращивать размер вектора

tic
log_of_random_vector = [];
for index = 1:length(random_vector)
    result_i = log(random_vector(index));
    log_of_random_vector = [log_of_random_vector, result_i];
end
disp(['duration of the very bad method on '...
    num2str(size_of_vector) ' elements'])
toc

Результат работы такой:

duration of the very bad method on 10000 elements Elapsed time is 0.286904 seconds.

Способ второй - заметно луче - заранее зарезервировать память под результат применения функции логарифма. Для наглядности увеличим размер вектора еще в 10 раз.

size_of_vector = 1000000;
random_vector = rand(1, size_of_vector);
 
log_of_random_vector = zeros(1, length(random_vector));
tic
for index = 1:length(random_vector)
    log_of_random_vector(index) = log(random_vector(index));
end
disp(['duration time of a bit improved method on '...
    num2str(size_of_vector) ' elements'])
toc

Результаты работы этого кода:

duration time of a bit improved method on 1000000 elements Elapsed time is 0.062692 seconds.

Способ третий - самый правильный и самый быстрый, работает мгновенно - взять функцию логарифма (то, что именно логарифма - не принципиально) сразу от вектора

tic
log_of_random_vector = log(random_vector);
disp(['duration of the best method on '...
    num2str(size_of_vector) ' elements'])
toc

Результат работы:

duration of the best method on 1000000 elements Elapsed time is 0.015231 seconds.

Пример 2.

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

То есть, если нам надо перемножить матрицы A и B, то не надо писать никаких циклов, а нужно просто написать C = A * B

A = randn(2000, 2000);
B = randn(2000, 2000);
tic
C = A * B;
disp('Matrix multiplication time')
toc

Результат: Matrix multiplication time Elapsed time is 0.612797 seconds.

Здесь надо заметить, что при использовании <<лобового>> метода умножения матриц размеров n * n время получения ответа равно O(n^3). Если предположить, что компьютер делает в секунду в среднем 10^6 операций, то что-то тут не сходится - получится, что компьютер за 0.612797 секунд выполнил 8*10^9 операции. Дело в том, что в MATLAB реализован быстрый алгоритм умножения матриц, в этом алгоритме сложность менее, чем кубическая от размера входных матриц. Поэтому при перемножении матриц важно пользоваться именно записью C = A * B, а не ударяться в написание тройного цикла - тот уж точно выполнит все 8*10^9 операций и хорошо, если за пару суток управится.

Пример 3. Есть еще несколько функций, которые заточены под матричные и векторные операции - это функции repmat и reshape. Остановимся на них подробнее.

Важно!!! Практика показала, что repmat и reshape в разных версиях MATLAB могут отличаться, поэтому будьте внимательны и осторожны при их использовании.

Пусть мы хотим сделать матрицу, в которой все строки одинаковые. Тут возможны несколько способов реализации. Способ первый - цикл.

random_row = randn(1, 1000);
n_rows = 100000;
repeating_row_matrix = zeros(n_rows, length(random_row));
tic
for index = 1:n_rows
    repeating_row_matrix(index, :) = random_row;
end
toc

Elapsed time is 1.596466 seconds.

Способ второй - использовать функцию repmat()

tic
repeating_row_matrix = repmat(random_row, n_rows, 1);
toc

Elapsed time is 0.311417 seconds.

Пусть теперь нам захотелось продублировать строки матрицы A, то есть вместо каждой строки матрицы написать две таких же строки, одну под другой. Как сделать это эффективно, то есть БЕЗ ЦИКЛОВ? Тут на помощь приходит функция reshape() (еще раз напоминаю про возможные различия в разных версиях).

A = [[1, 2, 3]; [4, 5, 6]; [7, 8, 9]];
disp(A)

Дальше ловкость рук и никакого мошенничества =)

B = reshape(repmat(A', 2, 1), 3, 6)';
disp(B)

Читателю предоставляется возможность разобраться с тем, как получилась матрица B, самостоятельно.

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

Пример 4, строим усеченную синусоиду, вариант 1

x = [0:0.01*2*pi:2*pi];	% the length of the vector is 101
y = sin(x);			% the 1st line
y1 = y;			% dummy for the 2nd line
CONST = 0.2;
for i = 1:length(y1)
	if y1(i) > CONST,	% please, do not act like that!
		y1(i) = CONST;
	elseif y1(i) < -CONST;
		y1(i) = -CONST;
	end
end
plot(x,y,'r-');		% plot the result
hold on
plot(x,y1,'b-');
hold off

Пример 4, вариант 2

n = 100;
x = linspace(0,2*pi,n);
y = sin(x);
idx = find(abs(y) > CONST);	% it's better to write
y(idx) = sign(y(idx))*CONST;	% here the same result in a few lines of code!

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

пример 5, вариант 1

n = 5;
x = rand(n,1);
A = magic(n);
for j = 1:n
	A(:,j) = A(:,j)*x(j);	% multiply some columns by some scalars
end

пример 5, вариант 2

n = 5;
x = rand(n,1);
A = magic(n);
A = A*diag(x);		% much better

Функции

функции бывают embedded, private, public, inline, см. doc; при написании желательно организовать код так, чтобы функция возвращала корректный результат независимо от того, что было подано в качестве аргумента — скаляр, вектор или матрица

пример, стандартизация вектора

function y = to01(x)
	y = (x - min(x))./(max(x) - min(x))
return
% place the code below in a separate file named 'to01.m'
% usage 1:
y = to01([1 2 3 4 5]);
% usage 2:
x = randn(1,20);
y = to01(x);
min(y)
max(y)

функции могут включать варьируемое число входных и выходных аргументов

function [a,b] = somefunc(c,d)
if nargin < 2, d = 1; end
% some code..., i.e.
a = c;
if nargout > 1, b = d; end
return

выполнение функций как текстовых строк в теле программы

eval('x = [1,2;3,4];')
[m,n] = eval('size(x)')
mat = feval('rand',2,2)

пример, как сделать случайную k-индексную матрицу

% how to make an [a x b x c x d x ... x e] (k-times)
k = 5;				% have to generate k-dimantional array
dims=randperm(k+1)		% let the dimension sizes be random
dims(find(dims==1)) = []	% remove the size of 1
str = sprintf('%d,',dims)	% type the vector into a string
str = sprintf('rand(%s);',str(1:end-1))	% add the rand function
mat = eval(str);		% evaluate the expression
size(mat)			% check the size of the obtained multi-dimensional array mat

инлайн-функции — функции задаваемые в теле программы

пример 1

f1 = inline('w(1) + w(2)*x','w','x');
y1 = f1([1 2],3)
y1 = f1([1,2],[1:5]')	% be care of this; use .-operations

пример 2

f2 = inline('(x - min(x))./(max(x) - min(x))','x');
y2 = f2([1:10])

Мини-тест

Мини-тест к Домашнему заданию-1 курса "Численные методы обучения по прецедентам"

• Шпаргалка

x = [1 3 5 7 9 11]'
x > 3 
find(x > 3)     % Find indices and values on non-zero elements.
any(x>3)        % Determine if any array elements are nonzero.
all(x>3)        % Determine if all array elements are nonzero.
sum(x)          % Sum of array elements.
cumsum(x)       % Cumulative sum.
diff(x)         % Differences
repmat(x,1,2)   % Replicate array.
reshape(x,2,3)	% Reshape array.
 
y = [11 13]'
intersect(x,y)	% Set intersection of two arrays.
ismember(x,y)	% Array elements that are members of set array.
issorted(x)     % Determine whether set elements are in sorted order.
setdiff(x,y)	% Set difference of two arrays.
setxor(x,y)     % Set exclusive OR of two arrays.
union(x,y)      % Set union of two arrays.
unique([x;y])	% Unique values in array

• Мини-тест, PDF
• Разбор теста, PDF

Рекомендации программистам

1. Matlab прост в освоении. Нет понятий проект, компилятор, библиотека.
2. Есть командная строка, редактор m-code, path list и help.
3. Если вы хотите использовать функцию из toolbox просто используйте ее имя в коде.
4. Коллизий имен не существует. Старое имя исчезает при его перезагрузке.
5. Функции вызываются по имени согласно path list (см. основное меню). Следите на очередностью в этом списке.
6. Почитайте help, если вам нужно написать известный алгоритм, возможно он уже написан.
7. В частности, функции, работающие с множествами: intersect, ismember, issorted, setdiff, setxor, union, unique и очень полезная функция is* находятся в Help navigator -> MATLAB -> Functions — Categorical list.
8. Wikipedia.org содержит много библиотек для Matlab.
9. Не нужно избегать сложных алгоритмов. Часто задачи классификации, регрессии, оптимизации решаются в одну строку кода.
10. Желательно документировать свои функции так:
    • первая строка — назначение функции
    • вторая строка — имя функции и ее входные и выходные аргументы
    • варианты использования функции
    • пример использования функции
11. Изучите все структуры данных в Matlab, их очень много. В последней версии (на момент написания основной части — R2008a) появилась поддержка работы с классами.
12. Избегайте циклов, если возможно, используйте операции работы с матрицами.
13. Имеется возможность подключать функции, написанные на C (в виде специально скомпилированных библиотек).

Некоммерческие версии

Matlab — коммерческая программа. Существуют некоммерческие варианты, совместимые по базовым конструкциям языка, но не совместимые по библиотечным функциям. Например, Scilab, Euler Math Toolbox и Octave.

События

Несколько раз в год фирма «Софтлайн» проводит семинары, посвященные новым версиям и отдельным подсистемам Matlab. Проводится Всероссийская научная конференция «Проектирование научных и инженерных приложений в среде MATLAB».

Смотри также

• MVR Composer — система порождения моделей нелинейной регрессии
• Автоматическая генерация документации
• Пример внешнего документирования систем
• Документирование функций Matlab
• Параллельные вычисления в Matlab

Внешние ссылки

• Matlab Webinars Лекции (голос и слайды, примерно пол-часа с возможностью перемотки) обычно состоят из введения в Matlab, решаемой задачи и рекламы. Рекомендуется к просмотру в том случае, когда нужно использовать какой-нибудь незнакомый toolbox. Можно увидеть интересный код или использование команд языка.
• Консультационный центр по Matlab
• Scilab, the open source platform for numerical computation
• Euler Math Toolbox
• Octave, high-level language for numerical computations
• Статья Matlab на Wikipedia.org, En. Смотри прочие ссылки по теме там.