Skip to content

Commit 23f05e5

Browse files
codearxivcodearxiv
authored
Add files via upload
1 parent d13d665 commit 23f05e5

18 files changed

+1933
-0
lines changed

‎MatchingPursuit.cpp

+78
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,78 @@
1+
// Copyright (C) 2019 Piotr (Peter) Beben <[email protected]>
2+
// See LICENSE included.
3+
4+
5+
#define EIGEN_NO_MALLOC
6+
7+
#include "MatchingPursuit.h"
8+
#include "ensure_buffer_size.h"
9+
#include "constants.h"
10+
11+
//-------------------------------------------------------------------------
12+
void MatchingPursuit::ensure(
13+
Index nd, Index na, Index lm)
14+
{
15+
if(nd != ndim || na != natm || lm != lmax){
16+
ndim = nd;
17+
natm = na;
18+
lmax = (na < lm) ? na : lm;
19+
}
20+
}
21+
22+
23+
//-------------------------------------------------------------------------
24+
/**
25+
Given a size n signal vector Y, an n x m dictionary matrix D of
26+
m atoms (unit column vectors of size n), and an integer L,
27+
our task is to find a sparse size m code vector X that encodes Y
28+
as closely as possible using no more than L atoms in D. In detail,
29+
we try to solve to minimization problem
30+
31+
min_X ||Y - DX||_F
32+
subject to
33+
||X||_0 <= L
34+
35+
where ||.||_F is the matrix Frobenius norm, and ||.||_0 is the
36+
vector L_0 norm (the number of non-zero entries in a vector).
37+
38+
This is a greedy approach using the matching pursuit algorithm.
39+
40+
@param[in] Y: size n vector.
41+
@param[in] D: n x m dictionary matrix.
42+
@param[in] latm: Sparsity constraint L.
43+
@param[out] X: size m code vector.
44+
@param[out] R: size n residual vector.
45+
46+
*/
47+
48+
49+
50+
void MatchingPursuit::operator()(
51+
const VectorXf& Y, const MatrixXf& D, Index latm, VectorXf& X, VectorXf& R)
52+
{
53+
assert(D.rows() == Y.rows());
54+
assert(D.cols() == X.rows());
55+
ensure(D.rows(), D.cols(), latm);
56+
57+
X.setZero();
58+
R = Y;
59+
60+
for(Index j = 1; j <= latm; ++j){
61+
float absprojmax = -float_infinity;
62+
float projmax = 0;
63+
Index imax = 0;
64+
for(Index i = 0; i < natm; ++i){
65+
if( X(i) != 0.0f ) continue;
66+
float proj = R.dot(D.col(i));
67+
float absproj = abs(proj);
68+
if( absproj > absprojmax ){
69+
projmax = proj;
70+
absprojmax = absproj;
71+
imax = i;
72+
}
73+
}
74+
X(imax) = projmax;
75+
R = R - projmax*D.col(imax);
76+
}
77+
}
78+
//-------------------------------------------------------------------------

‎MatchingPursuit.h

+39
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
//-----------------------------------------------------------
2+
// Copyright (C) 2019 Piotr (Peter) Beben <[email protected]>
3+
// See LICENSE included.
4+
5+
6+
#ifndef MATCHINGPURSUIT_H
7+
#define MATCHINGPURSUIT_H
8+
9+
#include "constants.h"
10+
11+
class MatchingPursuit
12+
{
13+
using Index = Eigen::Index;
14+
using MatrixXf = Eigen::MatrixXf;
15+
using VectorXf = Eigen::VectorXf;
16+
17+
public:
18+
MatchingPursuit() :
19+
ndim(0), natm(0), lmax(0) {}
20+
21+
MatchingPursuit(const MatchingPursuit &mp) : MatchingPursuit() {}
22+
//{ ensure(sa.ndim, sa.natm, sa.lmax); }
23+
24+
~MatchingPursuit(){}
25+
26+
void ensure(Index nd, Index na, Index lm);
27+
28+
void operator() (
29+
const VectorXf& Y, const MatrixXf& D, Index l,
30+
VectorXf& X, VectorXf& R);
31+
32+
private:
33+
Index ndim; // Size of signal vector.
34+
Index natm; // No. of 'atoms' in dictionary.
35+
Index lmax; // maximum sparsity constraint <= natm.
36+
37+
};
38+
39+
#endif // MATCHINGPURSUIT_H

‎OrthogonalPursuit.cpp

+173
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,173 @@
1+
// Copyright (C) 2019 Piotr (Peter) Beben <[email protected]>
2+
// See LICENSE included.
3+
4+
#define EIGEN_NO_MALLOC
5+
6+
#include "OrthogonalPursuit.h"
7+
#include "ensure_buffer_size.h"
8+
#include "constants.h"
9+
//#include <functional>
10+
11+
using Eigen::Matrix;
12+
using Eigen::Map;
13+
using Eigen::Dynamic;
14+
using Eigen::LDLT;
15+
using Eigen::Aligned16;
16+
17+
18+
19+
//-------------------------------------------------------------------------
20+
void OrthogonalPursuit::ensure(Index nd, Index na, Index lm)
21+
{
22+
if(nd != ndim || na != natm || lm != lmax){
23+
ndim = nd;
24+
natm = na;
25+
lmax = (na < lm) ? na : lm;
26+
ensureWorkspace();
27+
}
28+
}
29+
30+
//-------------------------------------------------------------------------
31+
32+
void OrthogonalPursuit::ensureWorkspace()
33+
{
34+
// Allocate more workspace as necessary.
35+
//std::function< size_t(Index) > align_padded =
36+
// [=](Index n) ->size_t { return ALIGNEDX*(1+(n/ALIGNEDX)); };
37+
size_t paddednd = align_padded(ndim);
38+
size_t paddedlm = align_padded(lmax);
39+
size_t paddedndlm = align_padded(ndim*lmax);
40+
size_t paddedlmsq = align_padded(lmax*lmax);
41+
42+
size_t ndworkNeed = paddednd + 3*paddedlm + 2*paddedndlm + 2*paddedlmsq;
43+
ensure_buffer_size(ndworkNeed+ndworkNeed/2, dwork);
44+
45+
// Map to pre-allocated workspace.
46+
size_t p = 0;
47+
new (&U) MapVectf(&dwork[0],ndim);
48+
p = p + paddednd;
49+
new (&V) MapVectf(&dwork[p],lmax);
50+
p = p + paddedlm;
51+
new (&W) MapVectf(&dwork[p],lmax);
52+
p = p + paddedlm;
53+
new (&XI) MapVectf(&dwork[p],lmax);
54+
p = p + paddedlm;
55+
new (&E) MapMtrxf(&dwork[p],ndim,lmax);
56+
p = p + paddedndlm;
57+
new (&F) MapMtrxf(&dwork[p],lmax,lmax);
58+
p = p + paddedlmsq;
59+
dworkOffset = p;
60+
61+
size_t niworkNeed = lmax;
62+
ensure_buffer_size(niworkNeed+niworkNeed/2, iwork);
63+
64+
// Map to pre-allocated workspace.
65+
new (&I) Map<Matrix<Index,Dynamic,1>>(&iwork[0],lmax);
66+
67+
if( lmax*lmax > ldltSize ){
68+
ldltSize = lmax*lmax;
69+
delete ldlt;
70+
ldlt = new LDLT<MatrixXf>(ldltSize);
71+
}
72+
}
73+
74+
//-------------------------------------------------------------------------
75+
/**
76+
Orthogonal Matching Pursuit:
77+
78+
Similar to matching pursuit, but provides a better approximation
79+
at the expense of significantly more computation. Namely, updates
80+
all the coefficients of the current code vector X' at each iteration
81+
so that DX' is an orthogonal projection of the signal vector Y onto
82+
the subspace spanned by the dictionary atoms corresponding to the
83+
nonzero entries of X'.
84+
85+
@param[in] Y: size n vector.
86+
@param[in] D: n x m dictionary matrix of unit column vectors.
87+
@param[in] latm: Sparsity constraint.
88+
@param[out] X: size m code vector.
89+
@param[out] R: size n residual vector.
90+
91+
*/
92+
93+
94+
95+
void OrthogonalPursuit::operator() (
96+
const VectorXf& Y, const MatrixXf& D, Index latm,
97+
VectorXf& X, VectorXf& R)
98+
{
99+
assert(D.rows() == Y.rows());
100+
assert(D.cols() == X.rows());
101+
ensure(D.rows(),D.cols(),latm);
102+
103+
X.setZero();
104+
R = Y;
105+
106+
for(Index j = 1; j <= latm; ++j){
107+
// Find the next 'nearest' atom to current residual R.
108+
float absprojmax = -float_infinity;
109+
float projmax = 0;
110+
Index imax = 0;
111+
for(Index i = 0; i < natm; ++i){
112+
if( X(i) != 0.0f ) continue;
113+
float proj = R.dot(D.col(i));
114+
float absproj = abs(proj);
115+
if( absproj > absprojmax ){
116+
projmax = proj;
117+
absprojmax = absproj;
118+
imax = i;
119+
}
120+
}
121+
U = D.col(imax); // Dictionary atom U 'nearest' to R
122+
E.col(j-1) = U; // ...save it in j^th column of E
123+
I(j-1) = imax; // ...and save column index of U
124+
X(imax) = 1.0f; // Set temporarily 1.0 to mark traversed.
125+
126+
// Map to pre-allocated workspace.
127+
Index p = dworkOffset;
128+
new (&ETblk) MapMtrxf(&dwork[p],j,ndim);
129+
p = p + align_padded(j*ndim);
130+
new (&Fblk) MapMtrxf(&dwork[p],j,j);
131+
132+
// With U added to the current set E of j nearest atoms,
133+
// optimise the coefficients of XI w.r.t this E. This is
134+
// done by projecting Y onto the subspace spanned by E.
135+
if( j > 1 ) {
136+
// Compute the product E^T(:,1:j) * E(:,1:j),
137+
// This can be done quicker by reusing the product
138+
// E^T(:,1:j-1) * E(:,1:j-1) from the previous
139+
// iteration.
140+
141+
V.segment(0,j-1).noalias() = U.transpose() * E.block(0,0,ndim,j-1);
142+
F.col(j-1).segment(0,j-1) = V.segment(0,j-1);
143+
F.row(j-1).segment(0,j-1) = V.segment(0,j-1);
144+
F(j-1,j-1) = 1.0f;
145+
146+
Fblk = F.block(0,0,j,j);
147+
ETblk = E.block(0,0,ndim,j).transpose();
148+
W.segment(0,j).noalias() = ETblk * Y;
149+
// Solve (E^T*E)*XI = (E^T)*Y
150+
ldlt->compute(Fblk);
151+
XI.segment(0,j) = ldlt->solve(W.segment(0,j));
152+
153+
//Update residual R
154+
R = Y;
155+
R.noalias() -= (E.block(0,0,ndim,j))*(XI.segment(0,j));
156+
}
157+
else{
158+
F(0,0) = 1.0f;
159+
XI(0) = Y.dot(U);
160+
//Update residual R
161+
R = Y;
162+
R.noalias() -= XI(0)*U;
163+
}
164+
}
165+
166+
// Map back to code vector.
167+
for(Index i = 0; i < latm; ++i){
168+
X(I(i)) = XI(i);
169+
}
170+
171+
172+
}
173+
//-------------------------------------------------------------------------

‎OrthogonalPursuit.h

+67
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,67 @@
1+
//-----------------------------------------------------------
2+
// Copyright (C) 2019 Piotr (Peter) Beben <[email protected]>
3+
// See LICENSE included.
4+
5+
#ifndef ORTHOGONALPURSUIT_H
6+
#define ORTHOGONALPURSUIT_H
7+
8+
9+
#include "constants.h"
10+
//#include <vector>
11+
12+
13+
class OrthogonalPursuit
14+
{
15+
using Index = Eigen::Index;
16+
using MatrixXf = Eigen::MatrixXf;
17+
using VectorXf = Eigen::VectorXf;
18+
using MapVectf = Eigen::Map<VectorXf, ALIGNEDX>;
19+
using MapMtrxf = Eigen::Map<MatrixXf, ALIGNEDX>;
20+
template<typename T> using Map = Eigen::Map<T, ALIGNEDX>;
21+
template<typename T> using vector = std::vector<T, Eigen::aligned_allocator<T>>;
22+
23+
public:
24+
OrthogonalPursuit() :
25+
ndim(0), natm(0), lmax(0), dworkOffset(0),
26+
U(nullptr,0), V(nullptr,0), W(nullptr,0),
27+
XI(nullptr,0), I(nullptr,0),
28+
E(nullptr,0,0), F(nullptr,0,0),
29+
ETblk(nullptr,0,0), Fblk(nullptr,0,0),
30+
ldlt(nullptr), ldltSize(0) {}
31+
32+
OrthogonalPursuit(const OrthogonalPursuit &op) : OrthogonalPursuit(){}
33+
//{ ensure(sa.ndim, sa.natm, sa.lmax); }
34+
35+
~OrthogonalPursuit(){
36+
delete ldlt;
37+
}
38+
39+
void ensure(Index nd, Index na, Index lm);
40+
41+
void operator() (
42+
const VectorXf& Y, const MatrixXf& D, Index l,
43+
VectorXf& X, VectorXf& R);
44+
45+
private:
46+
void ensureWorkspace();
47+
48+
Index ndim; // Size of signal vector.
49+
Index natm; // No. of 'atoms' in dictionary.
50+
Index lmax; // maximum sparsity constraint <= natm.
51+
52+
vector<float> dwork;
53+
size_t dworkOffset;
54+
vector<Index> iwork;
55+
// Work-space mapped onto work buffer.
56+
MapVectf U, V, W, XI;
57+
Map< Eigen::Matrix<Index,Eigen::Dynamic,1> > I;
58+
MapMtrxf E, F, ETblk, Fblk;
59+
60+
Eigen::LDLT<MatrixXf> *ldlt;
61+
Index ldltSize;
62+
63+
64+
};
65+
66+
#endif // ORTHOGONALPURSUIT_H
67+

‎constants.h

+26
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
//-----------------------------------------------------------
2+
// Copyright (C) 2019 Piotr (Peter) Beben <[email protected]>
3+
// See LICENSE included.
4+
5+
#ifndef CONSTANTS_H
6+
#define CONSTANTS_H
7+
8+
//#include <Eigen/Dense>
9+
//#include <limits>
10+
11+
#define __ALIGNED_MEMORY
12+
13+
#ifdef __ALIGNED_MEMORY
14+
const size_t ALIGNEDX = Eigen::Aligned16;
15+
inline size_t align_padded(size_t n) {
16+
return ALIGNEDX > 0 ? ALIGNEDX*(1+((n-1)/ALIGNEDX)) : n;
17+
}
18+
#else
19+
const size_t ALIGNEDX = Eigen::Unaligned;
20+
inline size_t align_padded(size_t n) { return n; }
21+
#endif
22+
23+
const float float_infinity = std::numeric_limits<float>::infinity();
24+
const double double_infinity = std::numeric_limits<double>::infinity();
25+
26+
#endif // CONSTANTS_H

0 commit comments

Comments
 (0)