1 // This is core/vnl/algo/vnl_svd.h
2 #ifndef vnl_svd_h_
3 #define vnl_svd_h_
4 //:
5 // \file
6 // \brief Holds the singular value decomposition of a vnl_matrix.
7 // \author Andrew W. Fitzgibbon, Oxford IERG
8 // \date 15 Jul 96
9 //
10 // \verbatim
11 // Modifications
12 // fsm, Oxford IESRG, 26 Mar 1999
13 // 1. The singular values are now stored as reals (not complexes) when T is complex.
14 // 2. Fixed bug : for complex T, matrices have to be conjugated as well as transposed.
15 // Feb.2002 - Peter Vanroose - brief doxygen comment placed on single line
16 // \endverbatim
17
18 #include <iosfwd>
19 #include <vnl/vnl_numeric_traits.h>
20 #include <vnl/vnl_vector.h>
21 #include <vnl/vnl_matrix.h>
22 #include <vnl/vnl_diag_matrix.h>
23 #include <vnl/algo/vnl_algo_export.h>
24 #ifdef _MSC_VER
25 # include <vcl_msvc_warnings.h>
26 #endif
27
28 //: Holds the singular value decomposition of a vnl_matrix.
29 //
30 // The class holds three matrices U, W, V such that the original matrix
31 // $M = U W V^\top$. The DiagMatrix W stores the singular values in decreasing
32 // order. The columns of U which correspond to the nonzero singular values
33 // form a basis for range of M, while the columns of V corresponding to the
34 // zero singular values are the nullspace.
35 //
36 // The SVD is computed at construction time, and inquiries may then be made
37 // of the SVD. In particular, this allows easy access to multiple
38 // right-hand-side solves without the bother of putting all the RHS's into a
39 // Matrix.
40 //
41 // This class is supplied even though there is an existing vnl_matrix method
42 // for several reasons:
43 //
44 // It is more convenient to use as it manages all the storage for
45 // the U,S,V matrices, allowing repeated queries of the same SVD
46 // results.
47 //
48 // It avoids namespace clutter in the Matrix class. While svd()
49 // is a perfectly reasonable method for a Matrix, there are many other
50 // decompositions that might be of interest, and adding them all would
51 // make for a very large Matrix class.
52 //
53 // It demonstrates the holder model of compute class, implementing an
54 // algorithm on an object without adding a member that may not be of
55 // general interest. A similar pattern can be used for other
56 // decompositions which are not defined as members of the library Matrix
57 // class.
58 //
59 // It extends readily to n-ary operations, such as generalized
60 // eigensystems, which cannot be members of just one matrix.
61
62 template <class T>
63 class VNL_ALGO_EXPORT vnl_svd
64 {
65 public:
66 //: The singular values of a matrix of complex<T> are of type T, not complex<T>
67 typedef typename vnl_numeric_traits<T>::abs_t singval_t;
68
69 //:
70 // Construct a vnl_svd<T> object from $m \times n$ matrix $M$. The
71 // vnl_svd<T> object contains matrices $U$, $W$, $V$ such that
72 // $U W V^\top = M$.
73 //
74 // Uses linpack routine DSVDC to calculate an ``economy-size'' SVD
75 // where the returned $U$ is the same size as $M$, while $W$
76 // and $V$ are both $n \times n$. This is efficient for
77 // large rectangular solves where $m > n$, typical in least squares.
78 //
79 // The optional argument zero_out_tol is used to mark the zero singular
80 // values: If nonnegative, any s.v. smaller than zero_out_tol in
81 // absolute value is set to zero. If zero_out_tol is negative, the
82 // zeroing is relative to |zero_out_tol| * sigma_max();
83
84 vnl_svd(vnl_matrix<T> const &M, double zero_out_tol = 0.0);
85 virtual ~vnl_svd() = default;
86
87 // Data Access---------------------------------------------------------------
88
89 //: find weights below threshold tol, zero them out, and update W_ and Winverse_
90 void zero_out_absolute(double tol = 1e-8); //sqrt(machine epsilon)
91
92 //: find weights below tol*max(w) and zero them out
93 void zero_out_relative(double tol = 1e-8); //sqrt(machine epsilon)
singularities()94 int singularities () const { return W_.rows() - rank(); }
rank()95 unsigned int rank () const { return rank_; }
well_condition()96 singval_t well_condition () const { return sigma_min()/sigma_max(); }
97
98 //: Calculate determinant as product of diagonals in W.
99 singval_t determinant_magnitude () const;
100 singval_t norm() const;
101
102 //: Return the matrix U.
U()103 vnl_matrix<T> & U() { return U_; }
104
105 //: Return the matrix U.
U()106 vnl_matrix<T> const& U() const { return U_; }
107
108 //: Return the matrix U's (i,j)th entry (to avoid svd.U()(i,j); ).
U(int i,int j)109 T U(int i, int j) const { return U_(i,j); }
110
111 //: Get at DiagMatrix (q.v.) of singular values, sorted from largest to smallest
W()112 vnl_diag_matrix<singval_t> & W() { return W_; }
113
114 //: Get at DiagMatrix (q.v.) of singular values, sorted from largest to smallest
W()115 vnl_diag_matrix<singval_t> const & W() const { return W_; }
Winverse()116 vnl_diag_matrix<singval_t> & Winverse() { return Winverse_; }
Winverse()117 vnl_diag_matrix<singval_t> const & Winverse() const { return Winverse_; }
W(int i,int j)118 singval_t & W(int i, int j) { return W_(i,j); }
W(int i)119 singval_t & W(int i) { return W_(i,i); }
sigma_max()120 singval_t sigma_max() const { return W_(0,0); } // largest
sigma_min()121 singval_t sigma_min() const { return W_(n_-1,n_-1); } // smallest
122
123 //: Return the matrix V.
V()124 vnl_matrix<T> & V() { return V_; }
125
126 //: Return the matrix V.
V()127 vnl_matrix<T> const& V() const { return V_; }
128
129 //: Return the matrix V's (i,j)th entry (to avoid svd.V()(i,j); ).
V(int i,int j)130 T V(int i, int j) const { return V_(i,j); }
131
132 //:
inverse()133 inline vnl_matrix<T> inverse () const { return pinverse(); }
134
135 //: pseudo-inverse (for non-square matrix) of desired rank.
136 vnl_matrix<T> pinverse (unsigned int rank = ~0u) const; // ~0u == (unsigned int)-1
137
138 //: Calculate inverse of transpose, using desired rank.
139 vnl_matrix<T> tinverse (unsigned int rank = ~0u) const; // ~0u == (unsigned int)-1
140
141 //: Recompose SVD to U*W*V', using desired rank.
142 vnl_matrix<T> recompose (unsigned int rank = ~0u) const; // ~0u == (unsigned int)-1
143
144 //: Solve the matrix equation M X = B, returning X
145 vnl_matrix<T> solve (vnl_matrix<T> const& B) const;
146
147 //: Solve the matrix-vector system M x = y, returning x.
148 vnl_vector<T> solve (vnl_vector<T> const& y) const;
149 void solve (T const *rhs, T *lhs) const; // min ||A*lhs - rhs||
150
151 //: Solve the matrix-vector system M x = y.
152 // Assuming that the singular values W have been preinverted by the caller.
153 void solve_preinverted(vnl_vector<T> const& rhs, vnl_vector<T>* out) const;
154
155 //: Return N such that M * N = 0
156 vnl_matrix<T> nullspace() const;
157
158 //: Return N such that M' * N = 0
159 vnl_matrix<T> left_nullspace() const;
160
161 //: Return N such that M * N = 0
162 vnl_matrix<T> nullspace(int required_nullspace_dimension) const;
163
164 //: Implementation to be done yet; currently returns left_nullspace(). - PVR.
165 vnl_matrix<T> left_nullspace(int required_nullspace_dimension) const;
166
167 //: Return the rightmost column of V.
168 // Does not check to see whether or not the matrix actually was rank-deficient -
169 // the caller is assumed to have examined W and decided that to his or her satisfaction.
170 vnl_vector<T> nullvector() const;
171
172 //: Return the rightmost column of U.
173 // Does not check to see whether or not the matrix actually was rank-deficient.
174 vnl_vector<T> left_nullvector() const;
175
valid()176 bool valid() const { return valid_; }
177
178 private:
179
180 int m_, n_; // Size of M, local cache.
181 vnl_matrix<T> U_; // Columns Ui are basis for range of M for Wi != 0
182 vnl_diag_matrix<singval_t> W_;// Singular values, sorted in decreasing order
183 vnl_diag_matrix<singval_t> Winverse_;
184 vnl_matrix<T> V_; // Columns Vi are basis for nullspace of M for Wi = 0
185 unsigned rank_;
186 bool have_max_;
187 singval_t max_;
188 bool have_min_;
189 singval_t min_;
190 double last_tol_;
191 bool valid_; // false if the NETLIB call failed.
192
193 // Disallow assignment.
vnl_svd(vnl_svd<T> const &)194 vnl_svd(vnl_svd<T> const &) { }
195 vnl_svd<T>& operator=(vnl_svd<T> const &) { return *this; }
196 };
197
198 template <class T>
199 inline
vnl_svd_inverse(vnl_matrix<T> const & m)200 vnl_matrix<T> vnl_svd_inverse(vnl_matrix<T> const& m)
201 {
202 return vnl_svd<T>(m).inverse();
203 }
204
205 template <class T>
206 std::ostream& operator<<(std::ostream&, vnl_svd<T> const& svd);
207
208 #endif // vnl_svd_h_
209