Re: Matrix inversion tool
On Sun, 9 Aug 2015, Oleg Goldshmidt wrote: Matan Ziv-Av ma...@svgalib.org writes: On Sat, 8 Aug 2015, Omer Zak wrote: What happens if you use the regular matrix inversion tool which works on real numbers? (after rejecting, as singular over Z2, matrices whose determinant modulo 2 is different from 1) This will work if the inverse only has integer entries. It may also work if your algorithm returns rational numbers, but will not work with floating point numbers (there are no reals in a computer). You cannot faithfully convert 1.1234 to an element of Z_2. The original matrix contains zeros and ones, and you never need to divide, just multiply and add modulo 2. Please read what you reply to. The suggestion was to treat the matrix as a real matrix, so the calculations will not be modulo 2. If you take the, for example, matrix [ [ 0, 0, 1, 1 ], [ 0, 1, 0, 1 ], [ 1, 0, 0, 1 ], [ 1, 1, 1, 0 ] ] The inverse, in floating point numbers, will be something like [ [ -0.1, -0.1, 0.3, 0.1 ], [ -0.1, 0.3, -0.1, 0.1 ], [ 0.3, -0.1, -0.1, 0.1 ], [ 0.1, 0.1, 0.1, -0.1 ] ] And you have no idea if 0.3 should be 0 or 1 mod 2. If you use rational numbers you get [ [ -1/3, -1/3, 2/3, 1/3 ], [ -1/3, 2/3, -1/3, 1/3 ], [ 2/3, -1/3, -1/3, 1/3 ], [ 1/3, 1/3, 1/3, -1/3 ] ] Which is easy to convert (you ignore the denominators, which will all be odd for an invertible matrix). -- Matan. ___ Linux-il mailing list Linux-il@cs.huji.ac.il http://mailman.cs.huji.ac.il/mailman/listinfo/linux-il
Re: Matrix inversion tool
On Sun, 9 Aug 2015, Oleg Goldshmidt wrote: Shachar Shemesh shac...@shemesh.biz writes: Hi all, I'm looking for a tool/code to invert a matrix. So far, this sounds trivial. I have one special requirement. I did not think it was too special, except I could not find anywhere that supplied it. I want the matrix to be over a different field (i.e. - not the real numbers). In my particular case, I want it to be over Z2 (zero and one). If you don't find ready code writing your own should be pretty simple for Z2 - fields don't get any simpler. ;-) In general, if you have a matrix A, then its inverse is given by A^(-1) = adj(A)/det(A) [1] The determinant of any matrix over Z2 is either 0 or 1. If it is 0 the matrix is not invertible. If it is 1, then the inverse matrix is the same as the adjugate of the original matrix. Now, the adjugate is the transpose of the cofactor matrix, which is C_ij = (-1)^(i+j)*M_ij[2] where M_ij is A's (ij)-minor, i.e., the determinant of the matrix derived by deleting the i-th row and the j-th column of A. Thus adj(A) = C_ji = (-1)^(i+j)*M_ji [3] Thus, from equations 1 and 3 we have {A^(-1)}_ij = (-1)^(i+j)*M_ji [4] There are probably library routines to compute minors and cofactors. If you don't find any, then, given that all you have is zeroes and ones, it should be easy and safe to write your own. Note that you need to compute the determinant to determine whether or not your matrix is invertible, and to compute the determinant you need the cofactors, since det(A) = Sum_{i=1}^{n} (A_ij*C_ij) [5] So, compute det(A), keep the cofactors, transpose (Eq. [3]) - done. Use The last line is wrong. The previous are mathematically correct, but computationally wrong. When we compute det(A) using [5], we only calculate the minors for the first line, so we only know the minors for the first line, we still need to calculate them for the rest of the lines. The calculation as suggested above will run in O(n!) steps, this will obviously stop working for matrices of size 15x15 or there about. The naive algorithm, taught to any engineering or science student in the first linear algebra course, is Gauss elimination (also known as LU decomposition in this context). It runs in O(n^3) steps. Note that this algorithm is used very similarly for both inverting a matrix and calculating the determinant. So even if we replace step [5] in your suggestion, with Gaussian elimination, we still get O(n^5), when our code actually has all necessary components for running in O(n^3) instead. This is why when teaching about the classical adjoint forumula for inverse ([1]), or equivalently, the Cramer's rule for solving linear equations, it is very important to impress upon the students that while it appears simpler than Gauss elimination, it is in fact a lot less efficient. As a side note, the best current algorithm runs in O(n^2.3727) (at least according to wikipedia), so if you care a lot about performance, you can still save some computation time by doing more research. -- Matan. ___ Linux-il mailing list Linux-il@cs.huji.ac.il http://mailman.cs.huji.ac.il/mailman/listinfo/linux-il
Re: Matrix inversion tool
Shachar Shemesh shac...@shemesh.biz writes: Hi all, I'm looking for a tool/code to invert a matrix. So far, this sounds trivial. I have one special requirement. I did not think it was too special, except I could not find anywhere that supplied it. I want the matrix to be over a different field (i.e. - not the real numbers). In my particular case, I want it to be over Z2 (zero and one). If you don't find ready code writing your own should be pretty simple for Z2 - fields don't get any simpler. ;-) In general, if you have a matrix A, then its inverse is given by A^(-1) = adj(A)/det(A) [1] The determinant of any matrix over Z2 is either 0 or 1. If it is 0 the matrix is not invertible. If it is 1, then the inverse matrix is the same as the adjugate of the original matrix. Now, the adjugate is the transpose of the cofactor matrix, which is C_ij = (-1)^(i+j)*M_ij[2] where M_ij is A's (ij)-minor, i.e., the determinant of the matrix derived by deleting the i-th row and the j-th column of A. Thus adj(A) = C_ji = (-1)^(i+j)*M_ji [3] Thus, from equations 1 and 3 we have {A^(-1)}_ij = (-1)^(i+j)*M_ji [4] There are probably library routines to compute minors and cofactors. If you don't find any, then, given that all you have is zeroes and ones, it should be easy and safe to write your own. Note that you need to compute the determinant to determine whether or not your matrix is invertible, and to compute the determinant you need the cofactors, since det(A) = Sum_{i=1}^{n} (A_ij*C_ij) [5] So, compute det(A), keep the cofactors, transpose (Eq. [3]) - done. Use normal integer arithmetic and take into account that over Z2 all odd numbers are 1 and all even numbers are 0 (i.e., get the integer result modulo 2). If you can overload + and * operators for type Z2 you can make it pretty simple. And efficient, if you implement addition as XOR and multiplication - as AND. If you find, e.g., a decent templatized C++ inversion routine maybe you can use it as is after overloading + and *? Hope it helps, -- Oleg Goldshmidt | p...@goldshmidt.org ___ Linux-il mailing list Linux-il@cs.huji.ac.il http://mailman.cs.huji.ac.il/mailman/listinfo/linux-il
Re: Matrix inversion tool
On Sun, 9 Aug 2015, Shachar Shemesh wrote: I should point out that the simplified matrix I was using to prove to myself that the idea has merit is 20x20 (inversed using state of the art in the vi text editing front). The code I'll actually have to run will be 273x273. The Matrix itself is sparsish, so it's really not that bad. Still, to prove things are working I'll have to (probably pre-calculate) several thousand such matrices. In fact, one of my problems in using tools such as Mathematica and friends is that I don't want to hand code the Matrix to be inversed. GAP (and similarly SAGE, maxima) can be scripted. You do not give your exact setup, but if you can write the matrices to a text file in any format, you can read it in GAP. If you can output to the same format as I posted, you do not even need to write any program in GAP to read it. -- Matan. ___ Linux-il mailing list Linux-il@cs.huji.ac.il http://mailman.cs.huji.ac.il/mailman/listinfo/linux-il
Re: Matrix inversion tool
On Sun, 9 Aug 2015, Shachar Shemesh wrote: I should point out that the simplified matrix I was using to prove to myself that the idea has merit is 20x20 (inversed using state of the art in the vi text editing front). The code I'll actually have to run will be 273x273. The Matrix itself is sparsish, so it's really not that bad. Still, to prove things are working I'll have to (probably pre-calculate) several thousand such matrices. Actually, for this specific problem, implementing Gaussian elimination is much better than looking for a library or a tool, and interfacing with it. Untested C-Pseudo code: Input: A,B: nxn matrices, B=I_n Output: det(A), B=A^{-1} if A det=1 for(i=0;in;i++) { for(j=i;jn;j++) if(A[j][i]==1) { ExchangeLines(A,i,j); ExchangeLines(B,i,j); break; } if(j==n) return 0; for(k=1;kn;k++) if(A[k][i]==1 k!=j) { AddLines(A,k,j); AddLines(B,k,j); } } return 1; -- Matan. ___ Linux-il mailing list Linux-il@cs.huji.ac.il http://mailman.cs.huji.ac.il/mailman/listinfo/linux-il
Re: Matrix inversion tool
On Aug 9, 2015 20:15, Matan Ziv-Av ma...@svgalib.org wrote: On Sun, 9 Aug 2015, Shachar Shemesh wrote: I should point out that the simplified matrix I was using to prove to myself that the idea has merit is 20x20 (inversed using state of the art in the vi text editing front). The code I'll actually have to run will be 273x273. The Matrix itself is sparsish, so it's really not that bad. Still, to prove things are working I'll have to (probably pre-calculate) several thousand such matrices. In fact, one of my problems in using tools such as Mathematica and friends is that I don't want to hand code the Matrix to be inversed. GAP (and similarly SAGE, maxima) can be scripted. You do not give your exact setup, but if you can write the matrices to a text file in any format, you can read it in GAP. If you can output to the same format as I posted, you do not even need to write any program in GAP to read it. Not being connected to the field of applied mathematics can only recommend solutions observed. There is vast amount of code already written for LinAlg. in python you can use numpy. In C(++) any BLAS library have all you need. OpenBLAS seems more effective. If yo need to manipulate matrices en masse strongly suggested to use GPU. Both OpenCL and CUDA are good. Best Regards, Evgeniy. ___ Linux-il mailing list Linux-il@cs.huji.ac.il http://mailman.cs.huji.ac.il/mailman/listinfo/linux-il