Movatterモバイル変換

[0]ホーム
URL:

Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly

 Sign up
Appearance settings

Commitbc4b4bc

Browse files
authored
Merge pull request#3169 from stan-dev/sum-to-zero-matrix
Add sum_to_zero_constrain and _free for matrix types
2 parents59c0014 +d49a254 commitbc4b4bc

File tree

5 files changed

+265
-70
lines changed

5 files changed

+265
-70
lines changed

‎stan/math/prim/constraint/sum_to_zero_constrain.hpp

Lines changed: 73 additions & 57 deletions
Original file line numberDiff line numberDiff line change
@@ -50,7 +50,7 @@ inline plain_type_t<Vec> sum_to_zero_constrain(const Vec& y) {
5050
value_type_t<Vec>sum_w(0);
5151
for (int i = N; i >0; --i) {
5252
double n =static_cast<double>(i);
53-
auto w =y_ref(i -1) *inv_sqrt(n * (n +1));
53+
auto w = y_ref.coeff(i -1) *inv_sqrt(n * (n +1));
5454
sum_w += w;
5555

5656
z.coeffRef(i -1) += sum_w;
@@ -61,47 +61,81 @@ inline plain_type_t<Vec> sum_to_zero_constrain(const Vec& y) {
6161
}
6262

6363
/**
64-
* Return avector with sumzerocorresponding to thespecified
65-
* freevector.
64+
* Return amatrix that sums tozeroover both therows
65+
*and columns corresponding to thefreematrix x.
6666
*
67-
* The sum-to-zero transform is defined using a modified version of the
68-
* the inverse of the isometric log ratio transform (ILR).
69-
* See:
70-
* Egozcue, Juan Jose; Pawlowsky-Glahn, Vera; Mateu-Figueras, Gloria;
71-
* Barcelo-Vidal, Carles (2003), "Isometric logratio transformations for
72-
* compositional data analysis", Mathematical Geology, 35 (3): 279–300,
73-
* doi:10.1023/A:1023818214614, S2CID 122844634
67+
* This is a linear transform, with no Jacobian.
7468
*
75-
* This implementation is closer to the description of the same using "pivot
76-
* coordinates" in
77-
* Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of
78-
* Compositional Data. In: Applied Compositional Data Analysis. Springer Series
79-
* in Statistics. Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3
69+
* @tparam Mat type of the matrix
70+
* @param x Free matrix input of dimensionality (N - 1, M - 1).
71+
* @return Zero-sum matrix of dimensionality (N, M).
72+
*/
73+
template<typename Mat,require_eigen_matrix_dynamic_t<Mat>* =nullptr,
74+
require_not_st_var<Mat>* =nullptr>
75+
inlineplain_type_t<Mat>sum_to_zero_constrain(const Mat& x) {
76+
constauto N = x.rows();
77+
constauto M = x.cols();
78+
79+
plain_type_t<Mat> Z =Eigen::MatrixXd::Zero(N +1, M +1);
80+
if (unlikely(N ==0 || M ==0)) {
81+
return Z;
82+
}
83+
auto&& x_ref =to_ref(x);
84+
85+
Eigen::Matrix<value_type_t<Mat>, -1,1> beta =Eigen::VectorXd::Zero(N);
86+
87+
for (int j = M -1; j >=0; --j) {
88+
value_type_t<Mat>ax_previous(0);
89+
90+
double a_j =inv_sqrt((j +1.0) * (j +2.0));
91+
double b_j = (j +1.0) * a_j;
92+
93+
for (int i = N -1; i >=0; --i) {
94+
double a_i =inv_sqrt((i +1.0) * (i +2.0));
95+
double b_i = (i +1.0) * a_i;
96+
97+
auto b_i_x = b_i * x_ref.coeff(i, j) - ax_previous;
98+
99+
Z.coeffRef(i, j) = (b_j * b_i_x) - beta.coeff(i);
100+
beta.coeffRef(i) += a_j * b_i_x;
101+
102+
Z.coeffRef(N, j) -= Z.coeff(i, j);
103+
Z.coeffRef(i, M) -= Z.coeff(i, j);
104+
105+
ax_previous += a_i * x_ref.coeff(i, j);
106+
}
107+
Z.coeffRef(N, M) -= Z.coeff(N, j);
108+
}
109+
110+
return Z;
111+
}
112+
113+
/**
114+
* Return a vector or matrix with sum zero corresponding to the specified
115+
* free input.
80116
*
81117
* This is a linear transform, with no Jacobian.
82118
*
83-
* @tparamVec type of the vector
119+
* @tparamT type of theinput, either avector or a matrix
84120
* @tparam Lp unused
85-
* @param y Free vectorinput of dimensionality K - 1.
121+
* @param y Free vectoror matrix
86122
* @param lp unused
87-
* @return Zero-sum vectorof dimensionality K.
123+
* @return Zero-sum vectoror matrix which is one larger in each dimension
88124
*/
89-
template<typename Vec,typename Lp,require_eigen_col_vector_t<Vec>* =nullptr,
90-
require_not_st_var<Vec>* =nullptr>
91-
inlineplain_type_t<Vec>sum_to_zero_constrain(const Vec& y, Lp& lp) {
125+
template<typename T,typename Lp, require_not_st_var<T>* =nullptr>
126+
inlineplain_type_t<T>sum_to_zero_constrain(const T& y, Lp& lp) {
92127
returnsum_to_zero_constrain(y);
93128
}
94129

95130
/**
96-
* Return a vector with sum zero corresponding to the specified
97-
* freevector.
131+
* Return a vectoror matrixwith sum zero corresponding to the specified
132+
* freeinput.
98133
* This overload handles looping over the elements of a standard vector.
99134
*
100-
* @tparam Vec A standard vector with inner type inheriting from
101-
* `Eigen::DenseBase` or a `var_value` with inner type inheriting from
102-
* `Eigen::DenseBase` with compile time dynamic rows and 1 column
103-
* @param[in] y free vector
104-
* @return Zero-sum vectors of dimensionality one greater than `y`
135+
* @tparam T A standard vector with inner type that is either a vector or a
136+
* matrix
137+
* @param[in] y free vector or matrix
138+
* @return Zero-sum vectors or matrices which are one larger in each dimension
105139
*/
106140
template<typename T,require_std_vector_t<T>* =nullptr>
107141
inlineautosum_to_zero_constrain(const T& y) {
@@ -114,13 +148,12 @@ inline auto sum_to_zero_constrain(const T& y) {
114148
* free vector.
115149
* This overload handles looping over the elements of a standard vector.
116150
*
117-
* @tparam Vec A standard vector with inner type inheriting from
118-
* `Eigen::DenseBase` or a `var_value` with inner type inheriting from
119-
* `Eigen::DenseBase` with compile time dynamic rows and 1 column
151+
* @tparam T A standard vector with inner type that is either a vector or a
152+
* matrix
120153
* @tparam Lp unused
121-
* @param[in] y free vector
154+
* @param[in] y free vector or matrix
122155
* @param[in, out] lp unused
123-
* @return Zero-sum vectorsof dimensionalityonegreater than `y`
156+
* @return Zero-sum vectorsor matrices which areonelarger in each dimension
124157
*/
125158
template<typename T,typename Lp,require_std_vector_t<T>* =nullptr,
126159
require_convertible_t<return_type_t<T>, Lp>* =nullptr>
@@ -130,36 +163,19 @@ inline auto sum_to_zero_constrain(const T& y, Lp& lp) {
130163
}
131164

132165
/**
133-
* Return a vector with sum zero corresponding to the specified
134-
* free vector.
135-
*
136-
* The sum-to-zero transform is defined using a modified version of
137-
* the inverse of the isometric log ratio transform (ILR).
138-
* See:
139-
* Egozcue, Juan Jose; Pawlowsky-Glahn, Vera; Mateu-Figueras, Gloria;
140-
* Barcelo-Vidal, Carles (2003), "Isometric logratio transformations for
141-
* compositional data analysis", Mathematical Geology, 35 (3): 279–300,
142-
* doi:10.1023/A:1023818214614, S2CID 122844634
143-
*
144-
* This implementation is closer to the description of the same using "pivot
145-
* coordinates" in
146-
* Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of
147-
* Compositional Data. In: Applied Compositional Data Analysis. Springer Series
148-
* in Statistics. Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3
149-
*
166+
* Return a vector or matrix with sum zero corresponding to the specified
167+
* free input.
150168
* This is a linear transform, with no Jacobian.
151169
*
152170
* @tparam Jacobian unused
153-
* @tparam Vec A type inheriting from `Eigen::DenseBase` or a `var_value` with
154-
* inner type inheriting from `Eigen::DenseBase` with compile time dynamic rows
155-
* and 1 column, or a standard vector thereof
171+
* @tparam T type of the input
156172
* @tparam Lp unused
157-
* @param[in] y free vector
173+
* @param[in] y free vector or matrix
158174
* @param[in, out] lp unused
159-
* @return Zero-sum vectorof dimensionalityonegreater than `y`
175+
* @return Zero-sum vectoror matrix which isonelarger in each dimension
160176
*/
161-
template<bool Jacobian,typenameVec,typename Lp>
162-
inlineplain_type_t<Vec>sum_to_zero_constrain(constVec& y, Lp& lp) {
177+
template<bool Jacobian,typenameT,typename Lp>
178+
inlineplain_type_t<T>sum_to_zero_constrain(constT& y, Lp& lp) {
163179
returnsum_to_zero_constrain(y);
164180
}
165181

‎stan/math/prim/constraint/sum_to_zero_free.hpp

Lines changed: 51 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55
#include<stan/math/prim/err.hpp>
66
#include<stan/math/prim/fun/Eigen.hpp>
77
#include<stan/math/prim/fun/to_ref.hpp>
8+
#include<stan/math/prim/fun/inv_sqrt.hpp>
89
#include<stan/math/prim/fun/sqrt.hpp>
910
#include<stan/math/prim/functor/apply_vector_unary.hpp>
1011
#include<cmath>
@@ -47,25 +48,70 @@ inline plain_type_t<Vec> sum_to_zero_free(const Vec& z) {
4748
return y;
4849
}
4950

50-
y.coeffRef(N -1) = -z_ref(N) *sqrt(N * (N +1)) / N;
51+
y.coeffRef(N -1) = -z_ref.coeff(N) *sqrt(N * (N +1)) / N;
5152

5253
value_type_t<Vec>sum_w(0);
5354

5455
for (int i = N -2; i >=0; --i) {
5556
double n =static_cast<double>(i +1);
56-
auto w =y(i +1) /sqrt((n +1) * (n +2));
57+
auto w = y.coeff(i +1) /sqrt((n +1) * (n +2));
5758
sum_w += w;
58-
y.coeffRef(i) = (sum_w -z_ref(i +1)) *sqrt(n * (n +1)) / n;
59+
y.coeffRef(i) = (sum_w - z_ref.coeff(i +1)) *sqrt(n * (n +1)) / n;
5960
}
6061

6162
return y;
6263
}
6364

6465
/**
65-
* Overload of `sum_to_zero_free()` to untransform each Eigen vector
66+
* Return an unconstrained matrix.
67+
*
68+
* @tparam Mat a column vector type
69+
* @param z Matrix of size (N, M)
70+
* @return Free matrix of length (N - 1, M - 1)
71+
* @throw std::domain_error if z does not sum to zero
72+
*/
73+
template<typename Mat,require_eigen_matrix_dynamic_t<Mat>* =nullptr>
74+
inlineplain_type_t<Mat>sum_to_zero_free(const Mat& z) {
75+
constauto& z_ref =to_ref(z);
76+
check_sum_to_zero("stan::math::sum_to_zero_free","sum_to_zero variable",
77+
z_ref);
78+
79+
constauto N = z_ref.rows() -1;
80+
constauto M = z_ref.cols() -1;
81+
82+
plain_type_t<Mat> x =Eigen::MatrixXd::Zero(N, M);
83+
if (unlikely(N ==0 || M ==0)) {
84+
return x;
85+
}
86+
87+
Eigen::Matrix<value_type_t<Mat>, -1,1> beta =Eigen::VectorXd::Zero(N);
88+
89+
for (int j = M -1; j >=0; --j) {
90+
value_type_t<Mat>ax_previous(0);
91+
92+
double a_j =inv_sqrt((j +1.0) * (j +2.0));
93+
double b_j = (j +1.0) * a_j;
94+
95+
for (int i = N -1; i >=0; --i) {
96+
double a_i =inv_sqrt((i +1.0) * (i +2.0));
97+
double b_i = (i +1.0) * a_i;
98+
99+
auto alpha_plus_beta = z_ref.coeff(i, j) + beta.coeff(i);
100+
101+
x.coeffRef(i, j) = (alpha_plus_beta + b_j * ax_previous) / (b_j * b_i);
102+
beta.coeffRef(i) += a_j * (b_i * x.coeff(i, j) - ax_previous);
103+
ax_previous += a_i * x.coeff(i, j);
104+
}
105+
}
106+
107+
return x;
108+
}
109+
110+
/**
111+
* Overload of `sum_to_zero_free()` to untransform each Eigen type
66112
* in a standard vector.
67113
* @tparam T A standard vector with with a `value_type` which inherits from
68-
* `Eigen::MatrixBase` with compile time rows or columns equal to 1.
114+
* `Eigen::MatrixBase`
69115
* @param z The standard vector to untransform.
70116
*/
71117
template<typename T,require_std_vector_t<T>* =nullptr>

‎stan/math/rev/constraint/sum_to_zero_constrain.hpp

Lines changed: 60 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55
#include<stan/math/rev/core/reverse_pass_callback.hpp>
66
#include<stan/math/rev/core/arena_matrix.hpp>
77
#include<stan/math/prim/fun/Eigen.hpp>
8+
#include<stan/math/prim/fun/inv_sqrt.hpp>
89
#include<stan/math/prim/fun/sqrt.hpp>
910
#include<stan/math/prim/constraint/sum_to_zero_constrain.hpp>
1011
#include<cmath>
@@ -55,14 +56,66 @@ inline auto sum_to_zero_constrain(T&& y) {
5556
double n =static_cast<double>(i +1);
5657

5758
// adjoint of the reverse cumulative sum computed in the forward mode
58-
sum_u_adj += arena_z.adj()(i);
59+
sum_u_adj += arena_z.adj().coeff(i);
5960

6061
// adjoint of the offset subtraction
61-
double v_adj = -arena_z.adj()(i +1) * n;
62+
double v_adj = -arena_z.adj().coeff(i +1) * n;
6263

6364
double w_adj = v_adj + sum_u_adj;
6465

65-
arena_y.adj()(i) += w_adj /sqrt(n * (n +1));
66+
arena_y.adj().coeffRef(i) += w_adj /sqrt(n * (n +1));
67+
}
68+
});
69+
70+
return arena_z;
71+
}
72+
73+
/**
74+
* Return a matrix that sums to zero over both the rows
75+
* and columns corresponding to the free matrix x.
76+
*
77+
* This is a linear transform, with no Jacobian.
78+
*
79+
* @tparam Mat type of the matrix
80+
* @param x Free matrix input of dimensionality (N - 1, M - 1).
81+
* @return Zero-sum matrix of dimensionality (N, M).
82+
*/
83+
template<typename T,require_rev_matrix_t<T>* =nullptr,
84+
require_not_t<is_rev_vector<T>>* =nullptr>
85+
inlineautosum_to_zero_constrain(T&& x) {
86+
using ret_type =plain_type_t<T>;
87+
if (unlikely(x.size() ==0)) {
88+
returnarena_t<ret_type>(Eigen::MatrixXd{{0}});
89+
}
90+
auto arena_x =to_arena(std::forward<T>(x));
91+
arena_t<ret_type> arena_z =sum_to_zero_constrain(arena_x.val());
92+
93+
reverse_pass_callback([arena_x, arena_z]()mutable {
94+
constauto Nf = arena_x.val().rows();
95+
constauto Mf = arena_x.val().cols();
96+
97+
Eigen::VectorXd d_beta =Eigen::VectorXd::Zero(Nf);
98+
99+
for (int j =0; j < Mf; ++j) {
100+
double a_j =inv_sqrt((j +1.0) * (j +2.0));
101+
double b_j = (j +1.0) * a_j;
102+
103+
double d_ax =0.0;
104+
105+
for (int i =0; i < Nf; ++i) {
106+
double a_i =inv_sqrt((i +1.0) * (i +2.0));
107+
double b_i = (i +1.0) * a_i;
108+
109+
double dY = arena_z.adj().coeff(i, j) - arena_z.adj().coeff(Nf, j)
110+
+ arena_z.adj().coeff(Nf, Mf) - arena_z.adj().coeff(i, Mf);
111+
double dI_from_beta = a_j * d_beta.coeff(i);
112+
d_beta.coeffRef(i) += -dY;
113+
114+
double dI_from_alpha = b_j * dY;
115+
double dI = dI_from_alpha + dI_from_beta;
116+
arena_x.adj().coeffRef(i, j) += b_i * dI + a_i * d_ax;
117+
d_ax -= dI;
118+
}
66119
}
67120
});
68121

@@ -89,12 +142,12 @@ inline auto sum_to_zero_constrain(T&& y) {
89142
*
90143
* This is a linear transform, with no Jacobian.
91144
*
92-
* @tparamVec type of the vector
93-
* @param y Free vectorinput of dimensionality K - 1.
145+
* @tparamT type of the vector or matrix
146+
* @param y Free vectoror matrix.
94147
* @param lp unused
95-
* @return Zero-sum vectorof dimensionality K.
148+
* @return Zero-sum vectoror matrix which is one larger in each dimension
96149
*/
97-
template<typename T,typename Lp,require_rev_col_vector_t<T>* =nullptr>
150+
template<typename T,typename Lp,is_rev_matrix<T>* =nullptr>
98151
inlineautosum_to_zero_constrain(T&& y, Lp& lp) {
99152
returnsum_to_zero_constrain(std::forward<T>(y));
100153
}

‎test/unit/math/mix/constraint/sum_to_zero_constrain_test.cpp

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -56,3 +56,24 @@ TEST(MathMixMatFun, sum_to_zeroTransform) {
5656
v5 <<1, -3,2,0, -1;
5757
sum_to_zero_constrain_test::expect_sum_to_zero_transform(v5);
5858
}
59+
60+
TEST(MathMixMatFun, sum_to_zero_matrixTransform) {
61+
Eigen::MatrixXdm0_0(0,0);
62+
sum_to_zero_constrain_test::expect_sum_to_zero_transform(m0_0);
63+
64+
Eigen::MatrixXdm1_1(1,1);
65+
m1_1 <<1;
66+
sum_to_zero_constrain_test::expect_sum_to_zero_transform(m1_1);
67+
68+
Eigen::MatrixXdm2_2(2,2);
69+
m2_2 <<1,2, -3,4;
70+
sum_to_zero_constrain_test::expect_sum_to_zero_transform(m2_2);
71+
72+
Eigen::MatrixXdm3_4(3,4);
73+
m3_4 <<1,2, -3,4,5,6, -7,8,9, -10,11, -12;
74+
75+
sum_to_zero_constrain_test::expect_sum_to_zero_transform(m3_4);
76+
77+
Eigen::MatrixXd m4_3 = m3_4.transpose();
78+
sum_to_zero_constrain_test::expect_sum_to_zero_transform(m4_3);
79+
}

0 commit comments

Comments
 (0)
[8]ページ先頭
©2009-2025 Movatter.jp