I tried the following to permute one matrix and vector:

Code: Select all

 
      typedef Eigen::Matrix<int, Eigen::Dynamic, 1, Eigen::ColMajor, DIM, 1> IndexType;

      template <typename Derived1, typename Derived2>
      void _solve_subproblem(Eigen::MatrixBase<Derived1> const & A, Eigen::MatrixBase<Derived2> const & b)
      {
         // m-th point has been added newly
         int m = A.rows();

         // Permute the data
         IndexType Z(m);            
         Z << 1, 2, 0, 3;

         // It seems that Eigen::PermutationMatrix does not give correct permuted matrix and vector.
         Eigen::PermutationMatrix<Eigen::Dynamic, DIM> P(Z);
         MatType Atilde = P*A*P.inverse();
         VecType btilde = P*b;


The expected result is something similar to:
Atilde = [ A(1, 1), A(1, 2), A(1,0), A(1, 3);
A(2, 1), A(2, 2), A(2,0), A(2, 3);
A(0, 1), A(0, 2), A(0,0), A(0, 3);
A(3, 1), A(3, 2), A(3, 0), A(3, 3)];
and
btilde = [b(1); b(2); b(0); b(3)].
But I got the wrong result, e.g.
btilde = [b(0); b(2); b(1); b(3)],
Atilde = [A(0, 0), A(0, 2), ......].

Do Eigen's permutation matrices work in the different way?
How can I do the intended operations?

Thanks in advance.

then you need to compute:

P.inverse()*A*P;
P.inverse()*b;

because P is the permutation:

Code: Select all

0 0 1 0
1 0 0 0
0 1 0 0
0 0 0 1


(each number indicates for each column the row number containing '1')

Thanks, ggael~

I did some experiments, and found out something worth commenting.

When I tested the following code:

Code: Select all

template <typename Derived1, typename Derived2>
void test_eigen_permutation(Eigen::MatrixBase<Derived1> const & A, Eigen::MatrixBase<Derived2> const & b)
{   
   std::cout << "A = [ " << A << " ]\n\n";
   std::cout << "b = [ " << b << " ]\n\n";

   Eigen::Vector4i Z;
   Z << 2, 1, 0, 3;
   //Z << 1, 2, 0, 3;
   std::cout << "Z = [ " << Z << " ]\n\n\n";


   // EIGEN FUNCTION
   Eigen::PermutationMatrix<Eigen::Dynamic, 4> P(Z);
   Eigen::Matrix4d Atilde = P*A*P.inverse();
   Eigen::Vector4d btilde = P*b;

   std::cout << "P*A*P^{-1} = [ " << Atilde << " ]\n\n";
   std::cout << "P*b = [ " << btilde << " ]\n\n\n";

   Atilde = P.inverse()*A*P;
   btilde = P.inverse()*b;

   std::cout << "P^{-1}*A*P = [ " << Atilde << " ]\n\n";
   std::cout << "P^{-1}*b = [ " << btilde << " ]\n\n\n";

   // STRAIGHT IMPLEMENTATION
   Eigen::Matrix4d A_;
   Eigen::Vector4d b_;
   for (int i = 0; i < 4; i++)
   {
      A_.row(i) = A.row(Z[i]);
      b_[i] = b[Z[i]];
   }

   for (int i = 0; i < 4; i++)
   {
      Atilde.col(i)  = A_.col(Z[i]);
      btilde[i] = b_[i];
   }

   std::cout << "permute(A) = [ " << Atilde << " ]\n\n";
   std::cout << "permute(b) = [ " << btilde << " ]\n\n\n";
}

int main(int argc, char* argv[])
{
   Eigen::Matrix4d A;
   A << 11, 12, 13, 14, 21, 22, 23, 24, 31, 32, 33, 34, 41, 42, 43, 44;

   Eigen::Vector4d b;
   b << 1, 2, 3, 4;

   test_eigen_permutation(A, b);

   return 0;
}


it seems that the result depends on the number of swaps of the permutation.
When the permutation includes one swap (e.g. Z << 2, 1, 0, 3;), I have the following result which shows that
everything works fine.

Code: Select all

A = [ 11 12 13 14
21 22 23 24
31 32 33 34
41 42 43 44 ]

b = [ 1
2
3
4 ]

Z = [ 2
1
0
3 ]


P*A*P^{-1} = [ 33 32 31 34
23 22 21 24
13 12 11 14
43 42 41 44 ]

P*b = [ 3
2
1
4 ]


P^{-1}*A*P = [ 33 32 31 34
23 22 21 24
13 12 11 14
43 42 41 44 ]

P^{-1}*b = [ 3
2
1
4 ]


permute(A) = [ 33 32 31 34
23 22 21 24
13 12 11 14
43 42 41 44 ]

permute(b) = [ 3
2
1
4 ]


However, when the permutation includes two swap (e.g. Z << 1, 2, 0, 3;), I have the following result which shows that
P*A*P.inverse() and P*b does not work fine.

A = [ 11 12 13 14
21 22 23 24
31 32 33 34
41 42 43 44 ]

b = [ 1
2
3
4 ]

Z = [ 1
2
0
3 ]


P*A*P^{-1} = [ 33 31 32 34
13 11 12 14
23 21 22 24
43 41 42 44 ]

P*b = [ 3
1
2
4 ]


P^{-1}*A*P = [ 22 23 21 24
32 33 31 34
12 13 11 14
42 43 41 44 ]

P^{-1}*b = [ 2
3
1
4 ]


permute(A) = [ 22 23 21 24
32 33 31 34
12 13 11 14
42 43 41 44 ]

permute(b) = [ 2
3
1
4 ]

ggael wrote:then you need to compute:

P.inverse()*A*P;
P.inverse()*b;

because P is the permutation:

Code: Select all

0 0 1 0
1 0 0 0
0 1 0 0
0 0 0 1


(each number indicates for each column the row number containing '1')

Another comments:

When the permutation includes three swaps, e.g. Z << 1, 2, 3, 0,
the result is as follows:

A = [ 11 12 13 14
21 22 23 24
31 32 33 34
41 42 43 44 ]

b = [ 1
2
3
4 ]

Z = [ 1
2
3
0 ]


P*A*P^{-1} = [ 44 41 42 43
14 11 12 13
24 21 22 23
34 31 32 33 ]

P*b = [ 4
1
2
3 ]


P^{-1}*A*P = [ 22 23 24 21
32 33 34 31
42 43 44 41
12 13 14 11 ]

P^{-1}*b = [ 2
3
4
1 ]


permute(A) = [ 22 23 24 21
32 33 34 31
42 43 44 41
12 13 14 11 ]

permute(b) = [ 2
3
4
1 ]

All these corroborates that you are correct. We have to use P*inverse()*A*P and P*inverse()*b always.

Thanks.

coolcat wrote:Thanks, ggael~

I did some experiments, and found out something worth commenting.

When I tested the following code:

Code: Select all

template <typename Derived1, typename Derived2>
void test_eigen_permutation(Eigen::MatrixBase<Derived1> const & A, Eigen::MatrixBase<Derived2> const & b)
{   
   std::cout << "A = [ " << A << " ]\n\n";
   std::cout << "b = [ " << b << " ]\n\n";

   Eigen::Vector4i Z;
   Z << 2, 1, 0, 3;
   //Z << 1, 2, 0, 3;
   std::cout << "Z = [ " << Z << " ]\n\n\n";


   // EIGEN FUNCTION
   Eigen::PermutationMatrix<Eigen::Dynamic, 4> P(Z);
   Eigen::Matrix4d Atilde = P*A*P.inverse();
   Eigen::Vector4d btilde = P*b;

   std::cout << "P*A*P^{-1} = [ " << Atilde << " ]\n\n";
   std::cout << "P*b = [ " << btilde << " ]\n\n\n";

   Atilde = P.inverse()*A*P;
   btilde = P.inverse()*b;

   std::cout << "P^{-1}*A*P = [ " << Atilde << " ]\n\n";
   std::cout << "P^{-1}*b = [ " << btilde << " ]\n\n\n";

   // STRAIGHT IMPLEMENTATION
   Eigen::Matrix4d A_;
   Eigen::Vector4d b_;
   for (int i = 0; i < 4; i++)
   {
      A_.row(i) = A.row(Z[i]);
      b_[i] = b[Z[i]];
   }

   for (int i = 0; i < 4; i++)
   {
      Atilde.col(i)  = A_.col(Z[i]);
      btilde[i] = b_[i];
   }

   std::cout << "permute(A) = [ " << Atilde << " ]\n\n";
   std::cout << "permute(b) = [ " << btilde << " ]\n\n\n";
}

int main(int argc, char* argv[])
{
   Eigen::Matrix4d A;
   A << 11, 12, 13, 14, 21, 22, 23, 24, 31, 32, 33, 34, 41, 42, 43, 44;

   Eigen::Vector4d b;
   b << 1, 2, 3, 4;

   test_eigen_permutation(A, b);

   return 0;
}


it seems that the result depends on the number of swaps of the permutation.
When the permutation includes one swap (e.g. Z << 2, 1, 0, 3;), I have the following result which shows that
everything works fine.

Code: Select all

A = [ 11 12 13 14
21 22 23 24
31 32 33 34
41 42 43 44 ]

b = [ 1
2
3
4 ]

Z = [ 2
1
0
3 ]


P*A*P^{-1} = [ 33 32 31 34
23 22 21 24
13 12 11 14
43 42 41 44 ]

P*b = [ 3
2
1
4 ]


P^{-1}*A*P = [ 33 32 31 34
23 22 21 24
13 12 11 14
43 42 41 44 ]

P^{-1}*b = [ 3
2
1
4 ]


permute(A) = [ 33 32 31 34
23 22 21 24
13 12 11 14
43 42 41 44 ]

permute(b) = [ 3
2
1
4 ]


However, when the permutation includes two swap (e.g. Z << 1, 2, 0, 3;), I have the following result which shows that
P*A*P.inverse() and P*b does not work fine.

A = [ 11 12 13 14
21 22 23 24
31 32 33 34
41 42 43 44 ]

b = [ 1
2
3
4 ]

Z = [ 1
2
0
3 ]


P*A*P^{-1} = [ 33 31 32 34
13 11 12 14
23 21 22 24
43 41 42 44 ]

P*b = [ 3
1
2
4 ]


P^{-1}*A*P = [ 22 23 21 24
32 33 31 34
12 13 11 14
42 43 41 44 ]

P^{-1}*b = [ 2
3
1
4 ]


permute(A) = [ 22 23 21 24
32 33 31 34
12 13 11 14
42 43 41 44 ]

permute(b) = [ 2
3
1
4 ]

ggael wrote:then you need to compute:

P.inverse()*A*P;
P.inverse()*b;

because P is the permutation:

Code: Select all

0 0 1 0
1 0 0 0
0 1 0 0
0 0 0 1


(each number indicates for each column the row number containing '1')

This is because in your special case "Z << 2, 1, 0, 3", the permutation is symmetric, and P = P^T = P^-1 .