Efficiency of matrix data operations

Specific Domains Data
1

Background:

I have the following matrix is raw data,for example:

a=\left[ \begin{matrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \\ &...& \end{matrix} \right] , b=\left[ \begin{matrix} 1 & 2 & 3 & 5 \\ 4 & 5 & 8 & 6 \\ 7 & 0 & 9 & 8 \\ 2 & 3 & 1 & 6 \\ 4 & 5 & 6 & 9 \\ 8 & 9 & 7 & 1 \\ &...&... \end{matrix} \right]

It can be discovered:(Also what I want to get through the program)

  • The first line of a : (\left[ \begin{matrix} 1 & 2 & 3 \end{matrix} \right]) is a subset of the first (\left[ \begin{matrix} 1 & 2 & 3 &5 \end{matrix} \right]) and fourth (\left[ \begin{matrix} 2 &3 &1&6 \end{matrix} \right]) lines of b.
  • The first line of b differs from the first line of a by “5”,
    the fourth line of b differs from the first line of a by “6”.
  • According to the above ideas, the result should be like:
c=\left[ \begin{matrix} 1 & 2 & 3 &1&4&5&6\\ 4 & 5 & 6 &2&5&8&9\\ 7 & 8 & 9&3&6&0&1\\ & & &...& & \end{matrix} \right]

My code:

c = hcat(a,zeros(Int64,size(a,1),4))    
    for i = 1:size(a,1)
            dex = Vector{Int64}()
            for j = 1:size(b,1)
                if a[i,:] ⊆ b[j,:]
                    dex = append!(dex,j)
                end
            end
            c[i,4] = dex[1]
            c[i,5] = dex[2]
    end 
    # Determine the number of lines a belongs to b,
    # here I am sure that there are only two lines in b containing a

    for i = 1:size(a,1)
            k1 = setdiff(b[a[i,4],:], a[i,:])
            k2 = setdiff(b[a[i,5],:], a[i,:])
            c[i,6] = k1[1]
            c[i,7] = k2[1]
    end
    # Write the element different from a in b to c

Question:

When my matrix size is large, the execution efficiency of the above code is really too low. :sweat_smile:
I want to ask if you have any good suggestions or other algorithms to achieve the same purpose.(That is to get the matrix c)
If you still need other information, please let me know in the comment, I will add it in time, thank you all. :smiley:

2

It is not clear to me what you are trying to do.

If you want to find repeated rows in a matrix or matrices, use an O(1) lookup mechanism, eg put the rows in a Dict.

3

I am actually surprise this does not throw error?
append!takes two collections,dexis a collection,j` is not?

edit: huh, it works on v1.0.3, I guess I am surprise I didn’t know that.
Btw, I don’t think append! need assignment

4

I think this is wrong, it’s really

        k1 = setdiff(b[c[i,4],:], a[i,:])
        k2 = setdiff(b[c[i,5],:], a[i,:])
5

I really suggest that you read the Performance Tips to improve your code.

I did not try to understand what you are trying to accomplish, but just by annotating the loops to avoid bound checking and using views, I could decrease the computational burden a lot:

function test1(a,b)
    c = hcat(a,zeros(Int64,size(a,1),4))
    for i = 1:size(a,1)
        dex = Vector{Int64}()
        for j = 1:size(b,1)
            if a[i,:] ⊆ b[j,:]
                dex = append!(dex,j)
            end
        end
        c[i,4] = dex[1]
        c[i,5] = dex[2]
    end
    # Determine the number of lines a belongs to b,
    # here I am sure that there are only two lines in b containing a

    for i = 1:size(a,1)
        k1 = setdiff(b[c[i,4],:], a[i,:])
        k2 = setdiff(b[c[i,5],:], b[i,:])
        c[i,6] = k1[1]
        c[i,7] = k2[1]
    end

    return c
end

function test2(a,b)
    c = hcat(a,zeros(Int64,size(a,1),4))

    @inbounds @views for i = 1:size(a,1)
        dex = Vector{Int64}()
        for j = 1:size(b,1)
            if a[i,:] ⊆ b[j,:]
                dex = push!(dex,j)
            end
        end
        c[i,4] = dex[1]
        c[i,5] = dex[2]
    end
    # Determine the number of lines a belongs to b,
    # here I am sure that there are only two lines in b containing a

    @inbounds for i = 1:size(a,1)
        k1 = setdiff(b[c[i,4],:], @view(a[i,:]))
        k2 = setdiff(b[c[i,5],:], @view(b[i,:]))
        c[i,6] = k1[1]
        c[i,7] = k2[1]
    end

    return c
end

P.S.: I modified the code as suggested by @cchderrick, otherwise it will not work.

Then I get:

julia> a = [1 2 3; 4 5 6; 7 8 9];

julia> b = [1 2 3 5;
            4 5 8 6;
            7 0 9 8;
            2 3 1 6;
            4 5 6 9;
            8 9 7 1];

julia> @btime test1($a,$b)
  3.490 ÎĽs (105 allocations: 9.83 KiB)
3Ă—7 Array{Int64,2}:
 1  2  3  1  4  5  6
 4  5  6  2  5  8  9
 7  8  9  3  6  0  1

julia> @btime test2($a,$b)
  2.318 ÎĽs (105 allocations: 7.20 KiB)
3Ă—7 Array{Int64,2}:
 1  2  3  1  4  5  6
 4  5  6  2  5  8  9
 7  8  9  3  6  0  1

and I am 100% sure that it can be improved. However, it will require knowledge of what the algorithm is doing.

2 Likes
6

yes! :sweat_smile:

7

If different iterations of your for loops are not related to each other use SIMD, to enhance the speed.
Check Matlab & Julia Matrix Operations Benchmark and particularly JuliaBenchSIMD for some inspiration on how to make use of SIMD. I have used SIMD in 3 MatrixAddition, MatrixMultiplication, ElementWiseOperations functions.

1 Like
8

Assuming that there are always two rows that match, no more and no less, and that size(a, 2) = 3 and size(b, 2) = 4, then I would replace the methods ⊆ and setdiff() by hard coded ones. You want to do this because if you do some “easy-profiling” like writing @time in front of the if a[i,:] ⊆ b[j,:] loops and in front of k1 = setdiff(b[c[i,4],:], a[i,:]) you will see that those are the lines where the script spends most of its time and almost all allocations. I also took the transpose the matrices a and b so that I’m accessing their elements column-wise (see the Performance Tips). @inbounds also helps.

My code is:

function foo2(a, b)
    c = zeros(Int, 4, size(a, 2))    
    @inbounds for i=1:size(a, 2)
        tmp = 1
        for j=1:size(b, 2)
            tmp > 2 && continue
            is_1 = a[1, i] == b[1, j] || a[2, i] == b[1, j] || a[3, i] == b[1, j]
            is_2 = a[1, i] == b[2, j] || a[2, i] == b[2, j] || a[3, i] == b[2, j]
            is_3 = a[1, i] == b[3, j] || a[2, i] == b[3, j] || a[3, i] == b[3, j]
            is_4 = a[1, i] == b[4, j] || a[2, i] == b[4, j] || a[3, i] == b[4, j]
            # This line could potentially save time with larger matrix sizes.
            #!is_1 & !is_2 & !is_3 & !is_4 && continue
            is_123 = is_1 & is_2 & is_3
            is_124 = is_1 & is_2 & is_4
            is_134 = is_1 & is_3 & is_4
            is_234 = is_2 & is_3 & is_4
            if is_123
                c[tmp, i] = j
                c[tmp+2, i] = b[4, j]
                tmp += 1
            elseif is_124
                c[tmp, i] = j
                c[tmp+2, i] = b[3, j]
                tmp += 1
            elseif is_134
                c[tmp, i] = j
                c[tmp+2, i] = b[2, j]
                tmp += 1
            elseif is_234
                c[tmp, i] = j
                c[tmp+2, i] = b[1, j]
                tmp += 1
            end
        end
    end
    return c
end

a = [1 2 3; 
     4 5 6;
     7 8 9];

b = [1 2 3 5;
     4 5 8 6;
     7 0 9 8;
     2 3 1 6;
     4 5 6 9;
     8 9 7 1];

aT = Matrix(transpose(a))
bT = Matrix(transpose(b))

Results:

>>> @btime foo($a, $b);
6.563 ÎĽs (103 allocations: 9.55 KiB)

>>> @btime foo2($a, $b);
173.424 ns (1 allocation: 176 bytes)

Nearly \times 40 faster. If you have problems with larger matrices then I’d try with multithreading in the outermost loop.

In any case, my script is not very elegant and is hard-coded for the input you gave, but maybe some of the above ideas help you.

1 Like
9

Maybe because I am not a native English speaker, there are some problems in expression, but your answer is what I want!!!:+1::+1::+1: Thank you very much :smile:.I hadn’t noticed before that “⊆” efficiency was very low.

Here I would like to continue to ask:

  • Are you talking about parallel processing at the outermost level like this: (because my actual matrix is really large)
using Distributed
addprocs(20)

function matrix_opt(a, b)
    c = zeros(Int, 4, size(a, 2))
    @inbounds @sync @distributed  for i=1:size(a, 2)
        tmp = 1
        for j=1:size(b, 2)
            tmp > 2 && continue
            is_1 = a[1, i] == b[1, j] || a[2, i] == b[1, j] || a[3, i] == b[1, j]
            is_2 = a[1, i] == b[2, j] || a[2, i] == b[2, j] || a[3, i] == b[2, j]
            is_3 = a[1, i] == b[3, j] || a[2, i] == b[3, j] || a[3, i] == b[3, j]
            is_4 = a[1, i] == b[4, j] || a[2, i] == b[4, j] || a[3, i] == b[4, j]
            # This line could potentially save time with larger matrix sizes.
            #!is_1 & !is_2 & !is_3 & !is_4 && continue
            is_123 = is_1 & is_2 & is_3
            is_124 = is_1 & is_2 & is_4
            is_134 = is_1 & is_3 & is_4
            is_234 = is_2 & is_3 & is_4
            if is_123
                c[tmp, i] = j
                c[tmp+2, i] = b[4, j]
                tmp += 1
            elseif is_124
                c[tmp, i] = j
                c[tmp+2, i] = b[3, j]
                tmp += 1
            elseif is_134
                c[tmp, i] = j
                c[tmp+2, i] = b[2, j]
                tmp += 1
            elseif is_234
                c[tmp, i] = j
                c[tmp+2, i] = b[1, j]
                tmp += 1
            end
        end
    end
    return Matrix(transpose(c))
end
  • I have another request similar to the title:
    I have the following two matrices:
a=\left[ \begin{matrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \\ &...& \end{matrix} \right] , b=\left[ \begin{matrix} 1&4&3&2\\ 7&2&9&8\\ 4&6&1&5\\ 8&2&1&7\\ 11&3&4&8\\ &...&... \end{matrix} \right]

I am quite sure:

  1. a are only three columns,b are only four columns.
  2. In b there must be only one row containing a.For example:
    the first line of a contains only the first line of b, and the second line of a contains only the third line of b.
  3. I need to return the line number contained in a and b .And I also need to get a line in b that contains any elements other than a.So what I end up with is that c should be like this:
c=\left[ \begin{matrix} 1 & 2 & 3 &1&4\\ 4 & 5 & 6 &3&1\\ 7 & 8 & 9 &2&2\\ &&...& \end{matrix} \right]

I would like to ask you to help me write a piece of code in accordance with this requirement.:blush:
And I hope it execution efficiency is high (consider multi-core parallelism).

Thank you very much. I can see your idea, but my reading ability in the code is not very strong. I hope to get your help. Thank you.:blush:

10

Hi, I think it was not your English, but the code that was a little hard to get. :slightly_smiling_face:

I would suggest @sync @distributed + SharedArrays or Threads.@threads in front of the for i=1:size(a, 2) loop. Please, try both and see which one works best for you.

I can see your idea, but my reading ability in the code is not very strong.

Many boolean expressions uncommented, sorry. The comments would be:

# True if the element b[2, j] is in a[:, i]. 
is_2 = a[1, i] == b[2, j] || a[2, i] == b[2, j] || a[3, i] == b[2, j]

...

# True if the elements of a[:, i] are contained in {b[1, j], b[3, j], b[4, j]}, no matter the order.
is_134 = is_1 & is_3 & is_4

...

# The element b[2, j] is not in a[:, i], so write b[2, j] in c.
elseif is_134
    c[tmp, i] = j
    c[tmp+2, i] = b[2, j]
    tmp += 1

I would tell you to try to make the code you asked for now that I explained the code a bit more. The code you want is very similar to the previous one.

In case you still don’t know how to make it here it is

using BenchmarkTools

function foo(a, b)
    c = zeros(Int, 2, size(a, 2))    
    @inbounds for i=1:size(a, 2)
        for j=1:size(b, 2)
            # True if the element b_{1, j} is in the ith row of a. 
            is_1 = a[1, i] == b[1, j] || a[2, i] == b[1, j] || a[3, i] == b[1, j]
            # True if the element b_{2, j} is in the ith row of a. 
            is_2 = a[1, i] == b[2, j] || a[2, i] == b[2, j] || a[3, i] == b[2, j]
            # True if ... 
            is_3 = a[1, i] == b[3, j] || a[2, i] == b[3, j] || a[3, i] == b[3, j]
            is_4 = a[1, i] == b[4, j] || a[2, i] == b[4, j] || a[3, i] == b[4, j]
            # This line could potentially save time with larger matrix sizes.
            #!is_1 & !is_2 & !is_3 & !is_4 && continue
            # True if the ith row of a is distributed in the first three elements of b[:, j]
            is_123 = is_1 & is_2 & is_3
            # True if the ith row of a is distributed in the elements: b[1, j], b[2, j], b[4, j].
            is_124 = is_1 & is_2 & is_4
            # True if the ith row of a is distributed in the elements: b[1, j], b[3, j], b[4, j].
            is_134 = is_1 & is_3 & is_4
            # True if ...
            is_234 = is_2 & is_3 & is_4
            if is_123
                c[1, i] = j
                c[2, i] = b[4, j]
            elseif is_124
                c[1, i] = j
                c[2, i] = b[3, j]
            elseif is_134
                c[1, i] = j
                c[2, i] = b[2, j]
            elseif is_234
                c[1, i] = j
                c[2, i] = b[1, j]
            end
        end
    end
    return c
end

a = [1 2 3; 4 5 6; 7 8 9];
b = [1 4 3 2; 7 2 9 8; 4 6 1 5; 8 2 1 7; 11 3 4 8];
aT = Matrix(transpose(a))
bT = Matrix(transpose(b))

c = foo(aT, bT)
display(hcat(a, c'))

@btime foo($aT, $bT);
11

I really appreciate your patience in answering this question in detail.:+1: Now I understand your code and my goal has been achieved. Thank you very much.:blush:

1 Like