多次元信号処理・テンソル計算用のMATLAB/Octave関数

MATLAB Octave 研究

ちまたでは，多次元配列信号処理と呼ばれるものがひっそりと研究されている．
多次元配列処理を行うためのMATLAB/Octave用の基礎的な関数（unfoldとfold，Nモード積）をつくったので，需要があるか分からないが晒しておく．

ダウンロードはこちらから．
http://dl.dropbox.com/u/525570/d.hatena/tensor_funcs.zip

多次元配列信号処理

多次元配列信号処理が普通の信号処理と異なる点は，多次元配列を扱うことである．
単一のマイクでとった音声はベクトルとして表すことができる．
複数のマイクでとった音声だったり，グレイスケール画象は行列で表すことができる．
ベクトルと行列と言えば線形代数ですよね．
したがって，線形代数の授業ででてきた行列式・一次従属・固有値分解・特異値分解といった手法に基づいて，信号を解析・変換・分解などすることができる．


多次元配列とは行列よりインデックスを持つ配列のことを指す．
例えば，ベクトル\b{x}の各要素はx_i，行列\b{X}の各要素はX_{ij}とかける．
この多次元配列では，要素を指す添字を増やすことができて，\mathcal{X}_{ijk\ldots}という感じになる．
ある要素を一つの立方体で表すと，ベクトル・行列・3次テンソルは以下のように描くことができる．

3次テンソルデータの例は動画（縦・横・時間という3つの配列軸を持ち，その中に画素値が入っている）などがある，
RGBの値をもつカラー画像も3次テンソルの中にデータを格納できる（縦・横・RGBという配列軸）．
4次以上のテンソルは，想像力豊かに想いを馳せて欲しい．
まぁ，4次テンソルの場合だと3次テンソルが何個かあるようなデータを想像してもらえると良いと思います．

1次，2次という数え方をすると，ベクトルの次元と混同しやすいよね．
英語だとN配列次元はN-wayとか言ったりする．
さらに，それぞれの配列次元の数え方は1モード，2モード，．．．という感じで呼ばれている．


細かいことを言い出すとボロが出そうなので，もっと知りたい方は，文献
T. G. Kolda, B. W. Bader, ``Tensor Decompositions and applications,'' SIAM review, vol. 51, no. 3, pp. 455-500, 2009.
を参照されたい．

Tensor toolbox

多次元配列処理をMATLABで行うために，Tensor toolboxというツールボックスがフリーで公開されている．

これで，テンソルを定義すると"Whos"をやったときの"Class"が"tensor"になる．
こんな感じ．

>> x = tenrand([10 20 30]);
>> whos
Name Size            Bytes    Class    Attributes
x        10x20x30  48376  tensor

さらに，いろいろなテンソル計算が定義されていて，テンソル分解のアルゴリズムもいくつか実装されている．

ただOctaveでは使えない．

テンソルをつくる

まず，テンソルをつくる．
実は，Tensor toolboxを導入しなくても最近のMATLAB/Octaveは多次元配列をサポートしている．
なので，普通に

x = randn([10 15 20]);
y = ones([10 15 20]);
z = zeros([10 15 20]);

とすれば，サイズ10×15×20の多次元配列（テンソル）ができる．

UnfoldingとFolding

多次元配列信号処理ではNモード積を多用する．
そのNモード積の計算に必要なのが，UnfoldingとFolding．
Unfoldingとはテンソルをあるモードに沿って並び替えて行列にする操作で，Foldingはその行列を再びテンソルに戻す．

Unfolding

unfold.m

function y = unfold(x, md)
% x: input tensor
% md: mode for unfolding
% y: output matrix

d = size(x);
if md == 1
	d1 = d(1);
	d2 = prod(d(2:end));
	y = reshape(x, d1, d2);
else
	rmode = setxor(1:length(d), md);
	z = permute(x, [md rmode]);
	y = reshape(z, d(md), prod(d(rmode)));
end

ifで場合分けしている理由は，モード1の展開の場合，permuteが要らないため．

Folding

fold.m

function y = fold(x, md, s)
% x: input matrix
% md: mode which the matrix is transposed along
% s: dimensions of output tensor
% y: output tensor

if md == 1 % if 1-mode is chosen for alonging
	y = reshape(x, s);
else
	Nd = length(s);
	rmode = setxor(1:length(s), md);
	z = reshape(x, [s(md) s(rmode)]);
	y = ipermute(z, [md rmode]);
end
テストプログラム

test2.m

clear all

%dim = randperm(3) + 1;
dim = [2 3 4];
fprintf('tensor size:');
fprintf(' %d x', dim);
fprintf('\b\b  \n');

x = reshape(randperm(prod(dim)), dim);
fprintf('tensor, x =\n')
disp(x)

% unfold and fold
for ii = 1:length(dim)
	fprintf('\n')
	fprintf('Unfolding along mode %d\n', ii)
	xn = unfold(x, ii);
	disp(xn)

	y = fold(xn, ii, dim);
	max_err = max(abs(vec(x - y)));
	fprintf('Error between original tensor and unfolded-folded tensor: %f\n', max_err)
end

Nモード積

Nモード積

mprod.m

function y = mprod(x, A, md);
% x: input tensor
% A: matrix for N-mode product
% md: mode
% y: N-mode producted tensor

Xn = unfold(x, md);
Yn = A'*Xn;

sz = size(x);
sz(md) = size(A, 2);

y = fold(Yn, md, sz);
テストプログラム

test3.m

clear all

dim = [2 3 4];
x = reshape(randperm(prod(dim)), dim);
fprintf('tensor (size:');
fprintf(' %d x', dim);
fprintf('\b\b):\n');

disp(x)

y = reshape(randperm(dim(2)*10), [dim(2) 10]);
fprintf('matrix (size:');
fprintf(' %d x %d):\n', dim(2), 10);
disp(y)

% mode product
z = mprod(x, y, 2);
dimt = size(z);
fprintf('2-mode product of x by y (size:');
fprintf(' %d x', dimt);
fprintf('\b\b):\n');
disp(z)

Tensor toolboxと実行速度の比較

Tensor toolboxと実行速度を比べるために以下のスクリプトを使った．
test_sp.m

clear all

addpath('tensor_toolbox')
addpath('tensor_toolbox/algorithms')
addpath('tensor_toolbox/met')

A = reshape(1:6, [3 2]);
dim = [2 3 4];
val = randperm(prod(dim));

%% tensor toolbox
x = tensor(val, dim);
y = ttm(x, A', 2);
fprintf('N mode product of tensor x and matrix A:\n')
disp(y)


clear x y

%% my functions and MATLAB built-in functions
x = reshape(val, dim);
y = mprod(x, A, 2);
fprintf('N mode product of tensor x and matrix A:\n')
disp(y)

clear all

%% calculation speed
by = 10.^(0:.2:2.6)
by = round(by);
Nr = 500;
t1 = zeros(length(by), Nr);
t2 = zeros(length(by), Nr);
u1 = zeros(length(by), Nr);
u2 = zeros(length(by), Nr);
for ii = 1:length(by)
	dim = [1 1 1]*by(ii);
	A = rand(dim(2), dim(2));

	for jj = 1:Nr;
		fprintf('%d\r', jj)

		tic;
		x = rand(dim);
		u1(ii, jj) = toc;
		tic
		y = mprod(x, A, 2);
		u2(ii, jj) = toc;

		clear x y

		tic;
		x = tenrand(dim);
		t1(ii, jj) = toc;
		tic
		y = ttm(x, A, 2);
		t2(ii, jj) = toc;

		clear x y
	end
	fprintf('tensor size:');
	fprintf(' %d x', dim);
	fprintf('\b\b  \n');
	fprintf(' - Creating random tensor:\n')
	fprintf('\t %1.3e \t %1.3e\n', mean(t1(ii, :)), mean(u1(ii, :)))
	fprintf(' - N-mode product:\n')
	fprintf('\t %1.3e \t %1.3e\n', mean(t2(ii, :)), mean(u2(ii, :)))
end

figure(1);
plot(log(by), log(mean(t1, 2)), '-ob', 'linewidth', 10, 'markersize', 10,...
	log(by), log(mean(u1, 2)), '-xr', 'linewidth', 10, 'markersize', 10)
legend('Tensor toolbox', 'My functions')
xlabel('dimensions, log_{10}(N)')
ylabel('calculation speed, log_{10}(time) [secs]')
print sp_rand.eps -depsc -F:20

figure(2);
plot(log(by), log(mean(t2, 2)), '-ob', 'linewidth', 10, 'markersize', 10,...
	log(by), log(mean(u2, 2)), '-xr', 'linewidth', 10, 'markersize', 10)
legend('Tensor toolbox', 'My functions')
xlabel('dimensions, log_{10}(N)')
ylabel('calculation speed, log_{10}(time) [secs]')
print sp_mprod.eps -depsc -F:20

ランダムテンソル（ランダムな値を持つテンソルを）を生成に要する時間と，テンソルと行列のNモード積に要する時間を測定した．
Tensor toolboxではttmという関数を用いた．
計算機はApple Xserve（Intel Xeon (2.93GHz)，12GB）で，MATLAB R2010a上で測定した．



これが，ランダムテンソル生成に要した時間



これが，Nモード積に要した時間

まぁ，似たようなもんですね．

その他