Brtdku

I'm trying to implement in C++ the calculation of the matrix permanent following the Glynn formula.

I try to explain briefly how this formula works. Suppose we have an nxn matrix.

| a b c |

| d e f | 

| g h i |

To calculate the permanent with the Glynn formula I should try to execute "a kind" of matrix product with a matrix that is a table of truths of length 2^n with n/2 rows and n columns.

Something like this. Suppose a matrix with n = 3.

| a b c | |+ + +|

| d e f | |+ - +|

| g h i | |+ + -|

          |+ - -|

development of the formula. I have to get:

∆(a + b + c)(d + e + f)(g + h + i)

where:

∆ è equal to the product of the first of the sign matrix (therefore + * + * + = +). The positive sign of a, b and c I got it by multiplying the first column of the matrix of letters by the sign of the first row of the sign matrix.

The order is then: multiply the first column of matrix A for the first row of the matrix of signs, multiply the second column of matrix A for the first row of the matrix of signs and then multiply the last column of matrix A for the first line of the matrix of the signs.

This is the first step. The second will be to multiply the first column of matrix A for the second row of the matrix of signs, multiply the second column of matrix A for the second row of the matrix of signs and then multiply the third column of matrix A for the second row of matrix signs and so on..

The end result is this:

= + (a + b + c)(d + e + f)(g + h + i) 

  - (a - b + c)(d - e + f)(g - h + i) 

  - (a + b - c)(d + e - f)(g + h - i)

  + (a - b - c)(d - e - f)(g - h - i)

I am trying to implement this algorithm in C. I have correctly created a random matrix nxn and the matrix of the signs I have encoded it in this way to perform the multiplications in a simple way (+ = 1 and - = -1).
So the product I have to do is:

| a b c | |1  1  1|

| d e f | |1 -1  1|

| g h i | |1  1 -1|

          |1 -1 -1|

The function that I tried to create to run those products and those sums is:

double permanent(double input_matrix[n][n], int sign_matrix[n]){



int rows =  pow (2, n) ;

int partial_result = 0;

int result = 1;



for(int r = 0; r < n; r++)

{

   for(int c = 0; c < n; c++)

   {

      partial = partial + input_matrix[c][r] * sign_matrix[r][c];

      //cout << parziale <<  endl;

   }



   cout << partial << endl;

   partial_result = partial_result * parziale;

   partial = 0;



}

The problem is that when I execute partial = partial + input_matrix [c][r] * sign_matrix [r][c]; I can not "hold the line of the matrix of signs steady, with the cause that the product is wrong because I multiply the first column of matrix A for the first row of the sign matrix (right). second row of the sign matrix (wrong! I should also multiply the second column for the first row of the sign matrix as in the written formula).

Any suggestions?

I think that you really need to check twice your understanding of the equation. Remember that each summation or multiplication is a for loop. So, by using the notation in the original paper, you get this:

#include <vector>

#include <cassert>



using mat = std::vector<std::vector<double>>;



double permanent(const mat& input_matrix, const mat& sign_matrix)

{

    int m = input_matrix.size();

    assert(m > 2);

    assert(input_matrix[0].size() == m);

    assert(sign_matrix.size() == m);

    int cols = sign_matrix[0].size();

    assert(cols == 1 << (m - 1));



    double result = 0;

    for (int t = 0; t < cols; ++t) {

        double delta = 1;

        for (int k = 0; k < m; ++k) {

            delta *= sign_matrix[k][t];

        }

        double p = 1;

        for (int j = 0; j < m; j++) {

            double s = 0;

            for (int i = 0; i < m; i++) {

                s += sign_matrix[i][t] * input_matrix[i][j];

            }

            p *= s;

        }

        result += delta * p;

    }

    return result / cols;

}



int main()

{

    mat A = { 

        { 1, 4, 7, }, 

        { 2, 5, 8, }, 

        { 3, 6, 9, },

    };

    mat sign_mat = {

        { 1,  1,  1,  1, },

        { 1, -1,  1, -1, },

        { 1,  1, -1, -1, },

    };

    auto perA = permanent(A, sign_mat);

}