Mentions légales du service

Skip to content
Snippets Groups Projects
Commit 1cda1295 authored by hhakim's avatar hhakim
Browse files

Fix matfaust for missing return of diagonal eigenvalues arg in fgft_palm (and... 

Fix matfaust for missing return of diagonal eigenvalues arg in fgft_palm (and mex function) and update a API doc.
parent 05a635df
Branches
Tags 2.4.1
No related merge requests found
Pipeline #833859 skipped
Hide whitespace changes
Inline Side-by-side
misc/test/src/Matlab/FaustFactoryTest.m
+ 210
− 210
View file @ 1cda1295
......@@ -18,162 +18,162 @@ classdef FaustFactoryTest < matlab.unittest.TestCase
end
methods (Test)
% function test_fact_palm4msa(this)
% disp('Test FaustFactory.fact_palm4msa()')
% import matfaust.*
% import matfaust.factparams.*
% num_facts = 2;
% is_update_way_R2L = false;
% init_lambda = 1.0;
% init_facts = cell(2,1);
% init_facts{1} = zeros(500,32);
% init_facts{2} = eye(32);
% %M = rand(500, 32)
% load([this.faust_paths{1},'../../../misc/data/mat/config_compared_palm2.mat']);
% % data matrix is loaded from file
% M = data;
% cons = cell(2,1);
% cons{1} = ConstraintInt(ConstraintName(ConstraintName.SPLIN), 500, 32, 5);
% cons{2} = ConstraintReal(ConstraintName(ConstraintName.NORMCOL), 32, 32, 1.0);
% stop_crit = StoppingCriterion(200);
% params = ParamsPalm4MSA(cons, stop_crit, 'is_update_way_R2L', is_update_way_R2L, 'init_lambda', init_lambda, 'step_size', ParamsFact.DEFAULT_STEP_SIZE,...
% 'constant_step_size', ParamsFact.DEFAULT_CONSTANT_STEP_SIZE);
% F = FaustFactory.fact_palm4msa(M, params)
% this.verifyEqual(size(F), size(M))
% %disp('norm F: ')
% %norm(F, 'fro')
% E = full(F)-M;
% disp('err: ')
% norm(E,'fro')/norm(M,'fro')
% % matrix to factorize and reference relative error come from
% % misc/test/src/C++/test_palm4MSA.cpp
% this.verifyEqual(norm(E,'fro')/norm(M,'fro'), 0.2710, 'AbsTol', 0.0001)
% end
%
% function test_fact_palm4msaCplx(this)
% disp('Test FaustFactory.fact_palm4msaCplx()')
% import matfaust.*
% import matfaust.factparams.*
% num_facts = 2;
% is_update_way_R2L = false;
% init_lambda = 1.0;
% init_facts = cell(2,1);
% init_facts{1} = zeros(500,32);
% init_facts{2} = eye(32);
% %M = rand(500, 32)
% load([this.faust_paths{1},'../../../misc/data/mat/config_compared_palm2.mat']);
% % data matrix is loaded from file
% M = data+j*data;
% cons = cell(2,1);
% cons{1} = ConstraintInt(ConstraintName(ConstraintName.SPLIN), 500, 32, 5);
% cons{2} = ConstraintReal(ConstraintName(ConstraintName.NORMCOL), 32, 32, 1.0);
% stop_crit = StoppingCriterion(200);
% params = ParamsPalm4MSA(cons, stop_crit, 'is_update_way_R2L', is_update_way_R2L, 'init_lambda', init_lambda);
% F = FaustFactory.fact_palm4msa(M, params)
% this.verifyEqual(size(F), size(M))
% %disp('norm F: ')
% %norm(F, 'fro')
% E = full(F)-M;
% disp('err: ')
% norm(E,'fro')/norm(M,'fro')
% this.verifyEqual(norm(E,'fro')/norm(M,'fro'), 0.29177, 'AbsTol', 0.0001)
% end
%
% function test_fact_hierarchical(this)
% disp('Test FaustFactory.fact_hierarchical()')
% import matfaust.*
% import matfaust.factparams.*
% num_facts = 4;
% is_update_way_R2L = false;
% init_lambda = 1.0;
% %init_facts = cell(num_facts,1);
% %init_facts{1} = zeros(500,32);
% %for i=2:num_facts
% %init_facts{i} = zeros(32);
% %end
% %M = rand(500, 32)
% load([this.faust_paths{1},'../../../misc/data/mat/matrix_hierarchical_fact.mat'])
% % matrix is loaded from file
% M = matrix;
% fact_cons = cell(3,1);
% res_cons = cell(3, 1);
% fact_cons{1} = ConstraintInt(ConstraintName(ConstraintName.SPLIN), 500, 32, 5);
% fact_cons{2} = ConstraintInt(ConstraintName(ConstraintName.SP), 32, 32, 96);
% fact_cons{3} = ConstraintInt(ConstraintName(ConstraintName.SP), 32, 32, 96);
% res_cons{1} = ConstraintReal(ConstraintName(ConstraintName.NORMCOL), 32, 32, 1);
% res_cons{2} = ConstraintInt(ConstraintName(ConstraintName.SP), 32, 32, 666);
% res_cons{3} = ConstraintInt(ConstraintName(ConstraintName.SP), 32, 32, 333);
% stop_crit = StoppingCriterion(200);
% stop_crit2 = StoppingCriterion(200);
% params = ParamsHierarchicalFact(fact_cons, res_cons, stop_crit, stop_crit2);
% F = FaustFactory.fact_hierarchical(M, params)
% this.verifyEqual(size(F), size(M))
% %disp('norm F: ')
% %norm(F, 'fro')
% E = full(F)-M;
% disp('err: ')
% norm(E,'fro')/norm(M,'fro')
% % matrix to factorize and reference relative error come from
% % misc/test/src/C++/hierarchicalFactorization.cpp
% this.verifyEqual(norm(E,'fro')/norm(M,'fro'), 0.26851, 'AbsTol', 0.00001)
% end
%
% function test_fact_hierarchicalCplx(this)
% disp('Test FaustFactory.fact_hierarchicalCplx()')
% import matfaust.*
% import matfaust.factparams.*
% num_facts = 4;
% is_update_way_R2L = false;
% init_lambda = 1.0;
% %init_facts = cell(num_facts,1);
% %init_facts{1} = zeros(500,32);
% %for i=2:num_facts
% %init_facts{i} = zeros(32);
% %end
% %M = rand(500, 32)
% load([this.faust_paths{1},'../../../misc/data/mat/matrix_hierarchical_fact.mat'])
% % matrix is loaded from file
% M = matrix+j*matrix;
% fact_cons = cell(3,1);
% res_cons = cell(3, 1);
% fact_cons{1} = ConstraintInt(ConstraintName(ConstraintName.SPLIN), 500, 32, 5);
% fact_cons{2} = ConstraintInt(ConstraintName(ConstraintName.SP), 32, 32, 96);
% fact_cons{3} = ConstraintInt(ConstraintName(ConstraintName.SP), 32, 32, 96);
% res_cons{1} = ConstraintReal(ConstraintName(ConstraintName.NORMCOL), 32, 32, 1);
% res_cons{2} = ConstraintInt(ConstraintName(ConstraintName.SP), 32, 32, 666);
% res_cons{3} = ConstraintInt(ConstraintName(ConstraintName.SP), 32, 32, 333);
% stop_crit = StoppingCriterion(200);
% stop_crit2 = StoppingCriterion(200);
% params = ParamsHierarchicalFact(fact_cons, res_cons, stop_crit, stop_crit2,...
% 'init_lambda', init_lambda, 'is_update_way_R2L', is_update_way_R2L);
% F = FaustFactory.fact_hierarchical(M, params)
% this.verifyEqual(size(F), size(M))
% %disp('norm F: ')
% %norm(F, 'fro')
% E = full(F)-M;
% disp('err: ')
% norm(E,'fro')/norm(M,'fro')
% % matrix to factorize and reference relative error come from
% % misc/test/src/C++/hierarchicalFactorization.cpp
% this.verifyEqual(norm(E,'fro')/norm(M,'fro'), 1.0063, 'AbsTol', 0.0001)
% end
function test_fact_palm4msa(this)
disp('Test FaustFactory.fact_palm4msa()')
import matfaust.*
import matfaust.factparams.*
num_facts = 2;
is_update_way_R2L = false;
init_lambda = 1.0;
init_facts = cell(2,1);
init_facts{1} = zeros(500,32);
init_facts{2} = eye(32);
%M = rand(500, 32)
load([this.faust_paths{1},'../../../misc/data/mat/config_compared_palm2.mat']);
% data matrix is loaded from file
M = data;
cons = cell(2,1);
cons{1} = ConstraintInt(ConstraintName(ConstraintName.SPLIN), 500, 32, 5);
cons{2} = ConstraintReal(ConstraintName(ConstraintName.NORMCOL), 32, 32, 1.0);
stop_crit = StoppingCriterion(200);
params = ParamsPalm4MSA(cons, stop_crit, 'is_update_way_R2L', is_update_way_R2L, 'init_lambda', init_lambda, 'step_size', ParamsFact.DEFAULT_STEP_SIZE,...
'constant_step_size', ParamsFact.DEFAULT_CONSTANT_STEP_SIZE);
F = FaustFactory.fact_palm4msa(M, params)
this.verifyEqual(size(F), size(M))
%disp('norm F: ')
%norm(F, 'fro')
E = full(F)-M;
disp('err: ')
norm(E,'fro')/norm(M,'fro')
% matrix to factorize and reference relative error come from
% misc/test/src/C++/test_palm4MSA.cpp
this.verifyEqual(norm(E,'fro')/norm(M,'fro'), 0.2710, 'AbsTol', 0.0001)
end
function test_fact_palm4msaCplx(this)
disp('Test FaustFactory.fact_palm4msaCplx()')
import matfaust.*
import matfaust.factparams.*
num_facts = 2;
is_update_way_R2L = false;
init_lambda = 1.0;
init_facts = cell(2,1);
init_facts{1} = zeros(500,32);
init_facts{2} = eye(32);
%M = rand(500, 32)
load([this.faust_paths{1},'../../../misc/data/mat/config_compared_palm2.mat']);
% data matrix is loaded from file
M = data+j*data;
cons = cell(2,1);
cons{1} = ConstraintInt(ConstraintName(ConstraintName.SPLIN), 500, 32, 5);
cons{2} = ConstraintReal(ConstraintName(ConstraintName.NORMCOL), 32, 32, 1.0);
stop_crit = StoppingCriterion(200);
params = ParamsPalm4MSA(cons, stop_crit, 'is_update_way_R2L', is_update_way_R2L, 'init_lambda', init_lambda);
F = FaustFactory.fact_palm4msa(M, params)
this.verifyEqual(size(F), size(M))
%disp('norm F: ')
%norm(F, 'fro')
E = full(F)-M;
disp('err: ')
norm(E,'fro')/norm(M,'fro')
this.verifyEqual(norm(E,'fro')/norm(M,'fro'), 0.29177, 'AbsTol', 0.0001)
end
function test_fact_hierarchical(this)
disp('Test FaustFactory.fact_hierarchical()')
import matfaust.*
import matfaust.factparams.*
num_facts = 4;
is_update_way_R2L = false;
init_lambda = 1.0;
%init_facts = cell(num_facts,1);
%init_facts{1} = zeros(500,32);
%for i=2:num_facts
%init_facts{i} = zeros(32);
%end
%M = rand(500, 32)
load([this.faust_paths{1},'../../../misc/data/mat/matrix_hierarchical_fact.mat'])
% matrix is loaded from file
M = matrix;
fact_cons = cell(3,1);
res_cons = cell(3, 1);
fact_cons{1} = ConstraintInt(ConstraintName(ConstraintName.SPLIN), 500, 32, 5);
fact_cons{2} = ConstraintInt(ConstraintName(ConstraintName.SP), 32, 32, 96);
fact_cons{3} = ConstraintInt(ConstraintName(ConstraintName.SP), 32, 32, 96);
res_cons{1} = ConstraintReal(ConstraintName(ConstraintName.NORMCOL), 32, 32, 1);
res_cons{2} = ConstraintInt(ConstraintName(ConstraintName.SP), 32, 32, 666);
res_cons{3} = ConstraintInt(ConstraintName(ConstraintName.SP), 32, 32, 333);
stop_crit = StoppingCriterion(200);
stop_crit2 = StoppingCriterion(200);
params = ParamsHierarchicalFact(fact_cons, res_cons, stop_crit, stop_crit2);
F = FaustFactory.fact_hierarchical(M, params)
this.verifyEqual(size(F), size(M))
%disp('norm F: ')
%norm(F, 'fro')
E = full(F)-M;
disp('err: ')
norm(E,'fro')/norm(M,'fro')
% matrix to factorize and reference relative error come from
% misc/test/src/C++/hierarchicalFactorization.cpp
this.verifyEqual(norm(E,'fro')/norm(M,'fro'), 0.26851, 'AbsTol', 0.00001)
end
function test_fact_hierarchicalCplx(this)
disp('Test FaustFactory.fact_hierarchicalCplx()')
import matfaust.*
import matfaust.factparams.*
num_facts = 4;
is_update_way_R2L = false;
init_lambda = 1.0;
%init_facts = cell(num_facts,1);
%init_facts{1} = zeros(500,32);
%for i=2:num_facts
%init_facts{i} = zeros(32);
%end
%M = rand(500, 32)
load([this.faust_paths{1},'../../../misc/data/mat/matrix_hierarchical_fact.mat'])
% matrix is loaded from file
M = matrix+j*matrix;
fact_cons = cell(3,1);
res_cons = cell(3, 1);
fact_cons{1} = ConstraintInt(ConstraintName(ConstraintName.SPLIN), 500, 32, 5);
fact_cons{2} = ConstraintInt(ConstraintName(ConstraintName.SP), 32, 32, 96);
fact_cons{3} = ConstraintInt(ConstraintName(ConstraintName.SP), 32, 32, 96);
res_cons{1} = ConstraintReal(ConstraintName(ConstraintName.NORMCOL), 32, 32, 1);
res_cons{2} = ConstraintInt(ConstraintName(ConstraintName.SP), 32, 32, 666);
res_cons{3} = ConstraintInt(ConstraintName(ConstraintName.SP), 32, 32, 333);
stop_crit = StoppingCriterion(200);
stop_crit2 = StoppingCriterion(200);
params = ParamsHierarchicalFact(fact_cons, res_cons, stop_crit, stop_crit2,...
'init_lambda', init_lambda, 'is_update_way_R2L', is_update_way_R2L);
F = FaustFactory.fact_hierarchical(M, params)
this.verifyEqual(size(F), size(M))
%disp('norm F: ')
%norm(F, 'fro')
E = full(F)-M;
disp('err: ')
norm(E,'fro')/norm(M,'fro')
% matrix to factorize and reference relative error come from
% misc/test/src/C++/hierarchicalFactorization.cpp
this.verifyEqual(norm(E,'fro')/norm(M,'fro'), 1.0063, 'AbsTol', 0.0001)
end
function test_fgft_givens(this)
disp('Test FaustFactory.fgft_givens()')
import matfaust.*
load([this.faust_paths{1} '../../../misc/data/mat/test_GivensDiag_Lap_U_J.mat'])
% Lap and J available
[F,D] = FaustFactory.fgft_givens(Lap, J)%, 0, 'verbosity', 1);
[F,D] = FaustFactory.fgft_givens(Lap, J);%, 0, 'verbosity', 1);
this.verifyEqual(size(F), size(Lap))
%disp('norm F: ')
%norm(F, 'fro')
E = F*diag(D)*F'-Lap;
E = F*full(D)*F'-Lap;
err = norm(E,'fro')/norm(Lap,'fro')
% the error reference is from the C++ test,
% misc/test/src/C++/GivensFGFT.cpp.in
this.verifyEqual(err, 0.0326529, 'AbsTol', 0.00001)
% verify it works the same using the eigtj() alias function
[F2,D2] = FaustFactory.eigtj(Lap, J)%, 0, 'verbosity', 2);
[F2,D2] = FaustFactory.eigtj(Lap, J);%, 0, 'verbosity', 2);
this.verifyEqual(full(F2),full(F))
this.verifyEqual(D,D2)
end
......@@ -184,86 +184,86 @@ classdef FaustFactoryTest < matlab.unittest.TestCase
load([this.faust_paths{1} '../../../misc/data/mat/test_GivensDiag_Lap_U_J.mat'])
% Lap and J available
t = size(Lap,1)/2;
[F,D] = FaustFactory.fgft_givens(Lap, J, t)%, 'verbosity', 2);
[F,D] = FaustFactory.fgft_givens(Lap, J, t); %, 'verbosity', 2);
this.verifyEqual(size(F), size(Lap))
%disp('norm F: ')
%norm(F, 'fro')
E = F*diag(D)*F'-Lap;
E = F*full(D)*F'-Lap;
err = norm(E,'fro')/norm(Lap,'fro')
% the error reference is from the C++ test,
% misc/test/src/C++/GivensFGFTParallel.cpp.in
this.verifyEqual(err,0.0410448, 'AbsTol', 0.00001)
% verify it works the same using the eigtj() alias function
[F2,D2] = FaustFactory.eigtj(Lap, J, t)%, 'verbosity', 2);
[F2,D2] = FaustFactory.eigtj(Lap, J, t); %, 'verbosity', 2);
this.verifyEqual(full(F2),full(F))
this.verifyEqual(D,D2)
end
% function test_fgft_palm(this)
% disp('Test FaustFactory.fgft_palm()')
% import matfaust.*
% import matfaust.factparams.*
% num_facts = 4;
% is_update_way_R2L = false;
% init_lambda = 1.0;
%
% load([this.faust_paths{1},'../../../misc/data/mat/HierarchicalFactFFT_test_U_L_params.mat'])
% % U, Lap, init_D, params are loaded from file
%
% fact_cons = cell(3,1);
% res_cons = cell(3, 1);
% fact_cons{1} = ConstraintInt(ConstraintName(ConstraintName.SP), 128, 128, 12288);
% fact_cons{2} = ConstraintInt(ConstraintName(ConstraintName.SP), 128, 128, 6144);
% fact_cons{3} = ConstraintInt(ConstraintName(ConstraintName.SP), 128, 128, 3072);
% res_cons{1} = ConstraintInt(ConstraintName(ConstraintName.SP), 128, 128, 384);
% res_cons{2} = ConstraintInt(ConstraintName(ConstraintName.SP), 128, 128, 384);
% res_cons{3} = ConstraintInt(ConstraintName(ConstraintName.SP), 128, 128, 384);
% stop_crit = StoppingCriterion(params.niter1);
% stop_crit2 = StoppingCriterion(params.niter2);
% params.fact_side = 0 % forced
% params.verbose = 0 % forced
% params.init_lambda = 128;
% params = ParamsHierarchicalFact(fact_cons, res_cons, stop_crit, stop_crit2, 'is_fact_side_left', params.fact_side == 1, 'is_update_way_R2L', params.update_way == 1, 'init_lambda', params.init_lambda, 'step_size', params.stepsize, 'constant_step_size', false, 'is_verbose', params.verbose ~= 1);
% diag_init_D = diag(init_D)
% [F,lambda] = FaustFactory.fgft_palm(U, Lap, params, diag_init_D)
% this.verifyEqual(size(F), size(U))
% %disp('norm F: ')
% %norm(F, 'fro')
% E = full(F)-U;
% err = norm(E,'fro')/norm(U,'fro')
% % matrix to factorize and reference relative error come from
% % misc/test/src/C++/hierarchicalFactorizationFFT.cpp
% this.verifyEqual(err, 0.05539, 'AbsTol', 0.00001)
% end
%
% function testHadamard(this)
% disp('Test FaustFactory.wht()')
% import matfaust.*
% n = randi([1,6]);
% H = FaustFactory.wht(n);
% fH = full(H);
% this.verifyEqual(nnz(fH), numel(fH));
% i = 1;
% while(i<2^n-1)
% j = i + 1;
% while(j<2^n)
% this.verifyEqual(fH(i,:)*fH(j,:)',0);
% j = j + 1;
% end
% i = i + 1;
% end
% end
%
% function testFourier(this)
% disp('Test FaustFactory.dft()')
% import matfaust.*
% n = randi([1,10]);
% F = FaustFactory.dft(n);
% fF = full(F);
% fftI = fft(eye(2^n));
% % this.verifyEqual(nnz(fH), numel(fH));
% this.verifyEqual(norm(fF-fftI), 0, 'AbsTol',10^-12);
% end
function test_fgft_palm(this)
disp('Test FaustFactory.fgft_palm()')
import matfaust.*
import matfaust.factparams.*
num_facts = 4;
is_update_way_R2L = false;
init_lambda = 1.0;
load([this.faust_paths{1},'../../../misc/data/mat/HierarchicalFactFFT_test_U_L_params.mat'])
% U, Lap, init_D, params are loaded from file
fact_cons = cell(3,1);
res_cons = cell(3, 1);
fact_cons{1} = ConstraintInt(ConstraintName(ConstraintName.SP), 128, 128, 12288);
fact_cons{2} = ConstraintInt(ConstraintName(ConstraintName.SP), 128, 128, 6144);
fact_cons{3} = ConstraintInt(ConstraintName(ConstraintName.SP), 128, 128, 3072);
res_cons{1} = ConstraintInt(ConstraintName(ConstraintName.SP), 128, 128, 384);
res_cons{2} = ConstraintInt(ConstraintName(ConstraintName.SP), 128, 128, 384);
res_cons{3} = ConstraintInt(ConstraintName(ConstraintName.SP), 128, 128, 384);
stop_crit = StoppingCriterion(params.niter1);
stop_crit2 = StoppingCriterion(params.niter2);
params.fact_side = 0 % forced
params.verbose = 0 % forced
params.init_lambda = 128;
params = ParamsHierarchicalFact(fact_cons, res_cons, stop_crit, stop_crit2, 'is_fact_side_left', params.fact_side == 1, 'is_update_way_R2L', params.update_way == 1, 'init_lambda', params.init_lambda, 'step_size', params.stepsize, 'constant_step_size', false, 'is_verbose', params.verbose ~= 1);
diag_init_D = diag(init_D)
[F,D,lambda] = FaustFactory.fgft_palm(U, Lap, params, diag_init_D)
this.verifyEqual(size(F), size(U))
%disp('norm F: ')
%norm(F, 'fro')
E = full(F)-U;
err = norm(E,'fro')/norm(U,'fro')
% matrix to factorize and reference relative error come from
% misc/test/src/C++/hierarchicalFactorizationFFT.cpp
this.verifyEqual(err, 0.05539, 'AbsTol', 0.00001)
end
function testHadamard(this)
disp('Test FaustFactory.wht()')
import matfaust.*
n = randi([1,6]);
H = FaustFactory.wht(n);
fH = full(H);
this.verifyEqual(nnz(fH), numel(fH));
i = 1;
while(i<2^n-1)
j = i + 1;
while(j<2^n)
this.verifyEqual(fH(i,:)*fH(j,:)',0);
j = j + 1;
end
i = i + 1;
end
end
function testFourier(this)
disp('Test FaustFactory.dft()')
import matfaust.*
n = randi([1,10]);
F = FaustFactory.dft(n);
fF = full(F);
fftI = fft(eye(2^n));
% this.verifyEqual(nnz(fH), numel(fH));
this.verifyEqual(norm(fF-fftI), 0, 'AbsTol',10^-12);
end
end
methods
......
wrapper/matlab/+matfaust/@FaustFactory/FaustFactory.m
+ 158
− 8
View file @ 1cda1295
......@@ -272,6 +272,82 @@ classdef FaustFactory
varargout = {F, lambda, p};
end
%===================================================================================
%> Computes the FGFT for the Fourier matrix U (it should be the eigenvectors of the Laplacian Lap).
%>
%> @note this algorithm is a variant of FaustFactory.fact_hierarchical.
%>
%> @param Lap The laplacian matrix.
%> @param U The Fourier matrix.
%> @param p The PALM hierarchical algorithm parameters.
%> @param init_D The initial diagonal vector. If none it will be the ones() vector by default.
%>
%> @retval [Uhat, Dhat, lambda, p]
%> - Uhat: the Faust factorization of U.
%> - Dhat: the diagonal matrix approximation of eigenvaules.
%> - lambda: see FaustFactory.fact_hierarchical
%> - p: see FaustFactory.fact_hierarchical
%>
%>
%> @b Example
%> @code
%> import matfaust.factparams.*
%>
%> % get the Laplacian
%> load('Laplacian_128_ring.mat');
%>
%> [U, D] = eig(Lap);
%> [D, I] = sort(diag(D));
%> D = diag(D);
%> U = U(:,I);
%>
%> dim = size(Lap, 1);
%>
%> nfacts = round(log2(dim)-3);
%> over_sp = 1.5; % sparsity overhead
%> dec_fact = .5; % decrease of the residuum sparsity
%>
%> % define the sparsity constraints for the factors
%> fact_cons = {};
%> res_cons = {};
%> for i=1:nfacts
%> fact_cons = [ fact_cons {ConstraintInt('sp', dim, dim, min(round(dec_fact^j*dim^2*over_sp), size(Lap,1)))} ];
%> res_cons = [ res_cons {ConstraintInt('sp', dim, dim, min(round(2*dim*over_sp), size(Lap, 1)))} ];
%> end
%>
%> % set the parameters for the PALM hierarchical algo.
%> params = ParamsHierarchicalFact(fact_cons, res_cons, StoppingCriterion(50), StoppingCriterion(100), 'step_size', 1e-6, 'constant_step_size', true, 'init_lambda', 1.0, 'is_fact_side_left', false);
%> %% compute FGFT for Lap, U, D
%> init_D_diag = diag(D);
%> [Uhat, Dhat, lambda, ~ ] = FaustFactory.fgft_palm(U, Lap, params, init_D_diag);
%>
%> %% errors on FGFT and Laplacian reconstruction
%> err_U = norm(Uhat-U, 'fro')/norm(U, 'fro')
%> err_Lap = norm(Uhat*full(Dhat)*Uhat'-Lap, 'fro') / norm(Lap, 'fro')
%> % Output:
%> % Faust::HierarchicalFact<FPP,DEVICE,FPP2>::compute_facts : factorization 1/4
%> % Faust::HierarchicalFact<FPP,DEVICE,FPP2>::compute_facts : factorization 2/4
%> % Faust::HierarchicalFact<FPP,DEVICE,FPP2>::compute_facts : factorization 3/4
%> % Faust::HierarchicalFact<FPP,DEVICE,FPP2>::compute_facts : factorization 4/4
%> % err_U =
%> % 1.0013
%> % err_Lap =
%> % 0.9707
%>
%> @endcode
%>
%> <p> @b See @b also FaustFactory.fact_hierarchical, FaustFactory.eigtj, FaustFactory.fgft_givens
%>
%> @b References:
%> - [1] Le Magoarou L., Gribonval R. and Tremblay N., "Approximate fast
%> graph Fourier transforms via multi-layer sparse approximations",
%> IEEE Transactions on Signal and Information Processing
%> over Networks 2018, 4(2), pp 407-420
%> <https://hal.inria.fr/hal-01416110>
%> - [2] Le Magoarou L. and Gribonval R., "Are there approximate Fast
%> Fourier Transforms on graphs ?", ICASSP, 2016. <https://hal.inria.fr/hal-01254108>
%>
%===================================================================================
function varargout = fgft_palm(U, Lap, p, varargin)
import matfaust.Faust
import matfaust.factparams.*
......@@ -331,18 +407,46 @@ classdef FaustFactory
% the setters for num_rows/cols verifies consistency with constraints
mex_params = struct('nfacts', p.num_facts, 'cons', {mex_constraints}, 'niter1', p.stop_crits{1}.num_its,'niter2', p.stop_crits{2}.num_its, 'sc_is_criterion_error', p.stop_crits{1}.is_criterion_error, 'sc_error_treshold', p.stop_crits{1}.error_treshold, 'sc_max_num_its', p.stop_crits{1}.max_num_its, 'sc_is_criterion_error2', p.stop_crits{2}.is_criterion_error, 'sc_error_treshold2', p.stop_crits{2}.error_treshold, 'sc_max_num_its2', p.stop_crits{2}.max_num_its, 'nrow', p.data_num_rows, 'ncol', p.data_num_cols, 'fact_side', p.is_fact_side_left, 'update_way', p.is_update_way_R2L, 'init_D', init_D, 'verbose', p.is_verbose, 'init_lambda', p.init_lambda);
if(isreal(U))
[lambda, core_obj] = mexHierarchical_factReal(U, mex_params, Lap);
[lambda, core_obj, Ddiag] = mexHierarchical_factReal(U, mex_params, Lap);
else
[lambda, core_obj] = mexHierarchical_factCplx(U, mex_params, Lap);
[lambda, core_obj, Ddiag] = mexHierarchical_factCplx(U, mex_params, Lap);
end
D = spdiag(Ddiag);
F = Faust(core_obj, isreal(U));
varargout = {F, lambda, p};
varargout = {F, D, lambda, p};
end
%==========================================================================================
%> @brief Computes the FGFT for the Laplacian matrix Lap.
%>
%>
%> @b Usage
%>
%> &nbsp;&nbsp;&nbsp; @b fgft_givens(Lap, J) calls the non-parallel Givens algorithm.<br/>
%> &nbsp;&nbsp;&nbsp; @b fgft_givens(Lap, J, 0) or fgft_givens(Lap, J, 1) do the same as in previous line.<br/>
%> &nbsp;&nbsp;&nbsp; @b fgft_givens(Lap, J, t) calls the parallel Givens algorithm (if t > 1, otherwise it calls basic Givens algorithm), see FaustFactory.eigtj. <br/>
%> &nbsp;&nbsp;&nbsp; @b fgft_givens(Lap, J, t, 'verbosity', 2) same as above with a level of verbosity of 2 in output. <br/>
%>
%> @param Lap the Laplacian matrix as a numpy array. Must be real and symmetric.
%> @param J see FaustFactory.eigtj
%> @param t see FaustFactory.eigtj
%> @param verbosity see FaustFactory.eigtj
%>
%> @retval [FGFT,D]:
%> - with FGFT being the Faust object representing the Fourier transform and,
%> - D as a sparse diagonal matrix of the eigenvalues in ascendant order along the rows/columns.
%>
%>
%> @b References:
%> - [1] Le Magoarou L., Gribonval R. and Tremblay N., "Approximate fast
%> graph Fourier transforms via multi-layer sparse approximations",
%> IEEE Transactions on Signal and Information Processing
%> over Networks 2018, 4(2), pp 407-420
%>
%>
%> <p> @b See @b also FaustFactory.eigtj, FaustFactory.fgft_palm
%>
%==========================================================================================
function [FGFT,D] = fgft_givens(Lap, J, varargin)
import matfaust.Faust
......@@ -354,12 +458,17 @@ classdef FaustFactory
if(size(Lap,1) ~= size(Lap,2))
error('Lap must be square')
end
if(~ isnumeric(J) || J-floor(J) > 0 || J <= 0)
error('J must be a positive integer.')
end
bad_arg_err = 'bad number of arguments.';
if(length(varargin) >= 1)
t = varargin{1};
if(~ isnumeric(t) || t-floor(t) > 0)
error('t must be an integer.')
if(~ isnumeric(t))
error('t must be a positive or nul integer.')
end
t = floor(abs(t));
t = min(t, J);
if(length(varargin) >= 2)
if(~ strcmp(varargin{2}, 'verbosity'))
error('arg. 4, if used, must be the str `verbosity''.')
......@@ -376,16 +485,57 @@ classdef FaustFactory
end
end
[core_obj, D] = mexfgftgivensReal(Lap, J, t, verbosity);
D = spdiag(D);
FGFT = Faust(core_obj, true);
end
%==========================================================================================
%> @brief Computes the eigenvalues and the eigenvectors transform (as a Faust object) using the truncated Jacobi algorithm.
%> @brief Computes the eigenvalues and the eigenvectors transform (as a Faust object) using the truncated Jacobi algorithm.
%>
%> The eigenvalues and the eigenvectors are approximate. The trade-off between accuracy and sparsity can be set through the parameters J and t.
%>
%>
%> @b Usage
%>
%> &nbsp;&nbsp;&nbsp; @b eigtj(M, J) calls the non-parallel Givens algorithm.<br/>
%> &nbsp;&nbsp;&nbsp; @b eigtj(M, J, 0) or eigtj(M, J, 1) do the same as in previous line.<br/>
%> &nbsp;&nbsp;&nbsp; @b eigtj(M, J, t) calls the parallel Givens algorithm (if t > 1, otherwise it calls basic Givens algorithm)<br/>
%> &nbsp;&nbsp;&nbsp; @b eigtj(M, J, t, 'verbosity', 2) same as above with a level of verbosity of 2 in output. <br/>
%>
%> @param M the matrix to diagonalize. Must be real and symmetric.
%> @param J defines the number of factors in the eigenvector transform V. The number of factors is round(J/t). Note that the last permutation factor is not in count here (in fact, the total number of factors in V is rather round(J/t)+1).
%> @param t the number of Givens rotations per factor. Note that t is forced to the value min(J,t). Besides, a value of t such that t > size(M,1)/2 won't lead to the desired effect because the maximum number of rotation matrices per factor is anyway size(M,1)/2. The parameter t is meaningful in the parallel version of the truncated Jacobi algorithm (cf. references below). If t <= 1 (by default) then the function runs the non-parallel algorithm.
%> @param verbosity the level of verbosity, the greater the value the more info. is displayed.
%>
%> @retval [V,W]
%> - V the Faust object representing the approximate eigenvector transform. V has its last factor being a permutation matrix, the goal of this factor is to apply to the columns of V the same order as eigenvalues set in W.
%> - W the approximate sparse diagonal matrix of the eigenvalues (in ascendant order along the rows/columns).
%>
%> @b Example
%> @code
%> import matfaust.*
%>
%> % get a Laplacian to diagonalize
%> load('Laplacian_256_community.mat')
%> % do it
%> [Uhat, Dhat] = FaustFactory.eigtj(Lap, size(Lap,1)*100, size(Lap, 1)/2, 'verbosity', 2)
%> % Uhat is the Fourier matrix/eigenvectors approximattion as a Faust (200 factors + permutation mat.)
%> % Dhat the eigenvalues diagonal matrix approx.
%> @endcode
%>
%>
%>
%> @b References:
%> - [1] Le Magoarou L., Gribonval R. and Tremblay N., "Approximate fast
%> graph Fourier transforms via multi-layer sparse approximations",
%> IEEE Transactions on Signal and Information Processing
%> over Networks 2018, 4(2), pp 407-420
%>
%> <p> @b See @b also FaustFactory.fgft_givens, FaustFactory.fgft_palm
%>
%==========================================================================================
function [FGFT,D] = eigtj(Lap, J, varargin)
[FGFT, D] = matfaust.FaustFactory.fgft_givens(Lap, J, varargin{:})
function [V,W] = eigtj(M, J, varargin)
[V, W] = matfaust.FaustFactory.fgft_givens(M, J, varargin{:});
end
%==========================================================================================
......
wrapper/matlab/src/mexHierarchical_fact.cpp.in
+ 7
− 0
View file @ 1cda1295
......@@ -138,6 +138,13 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
Faust::TransformHelper<SCALAR,Cpu>* F = new Faust::TransformHelper<SCALAR, Cpu>(t);
plhs[1] = convertPtr2Mat<Faust::TransformHelper<SCALAR, Cpu>>(F);
if(is_fgft)
{
Faust::Vect<SCALAR,Cpu> D(matrix.getNbRow());
(dynamic_cast<Faust::HierarchicalFactFFT<SCALAR,Cpu,FPP2>*>(hier_fact))->get_D(const_cast<SCALAR*>(D.getData())); // not respecting constness for optimiation (saving a vector copy)
plhs[2] = FaustVec2mxArray(D);
}
delete hier_fact;
}
catch (const std::exception& e)
......
wrapper/python/pyfaust.py
+ 11
− 9
View file @ 1cda1295
......@@ -1836,6 +1836,7 @@ class FaustFactory:
# errors on FGFT and Laplacian reconstruction
err_U = (Uhat-U).norm()/norm(U)
err_Lap = norm(Uhat.todense()*diag(Dhat)*Uhat.T.todense()-Lap)/norm(Lap)
err_Lap = norm(Uhat.todense()*diag(Dhat)*Uhat.T.todense())/norm(Lap)
print("err_U:", err_U)
print("err_Lap:", err_Lap)
......@@ -1846,7 +1847,7 @@ class FaustFactory:
# Faust::HierarchicalFact<FPP,DEVICE,FPP2>::compute_facts : factorization 2/3
# Faust::HierarchicalFact<FPP,DEVICE,FPP2>::compute_facts : factorization 3/3
# err_U: 1.00138959974
# err_Lap: 0.00319004898686
# err_Lap: 0.997230709335
Args:
......@@ -1874,8 +1875,8 @@ class FaustFactory:
References:
- [1] Le Magoarou L., Gribonval R. and Tremblay N., "Approximate fast
graph Fourier transforms via multi-layer sparse approximations",
submitted to IEEE Transactions on Signal and Information Processing
over Networks.
IEEE Transactions on Signal and Information Processing
over Networks 2018, 4(2), pp 407-420
<https://hal.inria.fr/hal-01416110>
- [2] Le Magoarou L. and Gribonval R., "Are there approximate Fast
Fourier Transforms on graphs ?", ICASSP, 2016. <https://hal.inria.fr/hal-01254108>
......@@ -1940,8 +1941,8 @@ class FaustFactory:
References:
[1] Le Magoarou L., Gribonval R. and Tremblay N., "Approximate fast
graph Fourier transforms via multi-layer sparse approximations",
submitted to IEEE Transactions on Signal and Information Processing
over Networks.
IEEE Transactions on Signal and Information Processing
over Networks 2018, 4(2), pp 407-420
<https://hal.inria.fr/hal-01416110>
Example:
......@@ -1957,8 +1958,9 @@ class FaustFactory:
Lap = data_dict['Lap'].astype(np.float)
Uhat, Dhat = FF.eigtj(Lap,J=Lap.shape[0]*100,t=int(Lap.shape[0]/2))
# Uhat is the Faust Fourier matrix approximation (200 factors + permutation mat.)
# Dhat the eigenvalues diagonal approx.
# Uhat is the Fourier matrix/eigenvectors approximattion as a Faust
# (200 factors + permutation mat.)
# Dhat the eigenvalues diagonal matrix approx.
print(Uhat)
print(Dhat)
......@@ -1989,8 +1991,8 @@ class FaustFactory:
References:
[1] Le Magoarou L., Gribonval R. and Tremblay N., "Approximate fast
graph Fourier transforms via multi-layer sparse approximations",
submitted to IEEE Transactions on Signal and Information Processing
over Networks.
IEEE Transactions on Signal and Information Processing
over Networks 2018, 4(2), pp 407-420.
<https://hal.inria.fr/hal-01416110>
<b/> See also FaustFactory.eigtj, FaustFactory.fgft_palm
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment