Embarrassingly, there is yet another bug.
while sp.colptr[col+1]<i
should be
while sp.colptr[col+1]<=i
(Arrg... where is the unsend/edit button)
In case someone ever, copy-pastes this. The (currently) correct version is:
function AtrA(sp::SparseMatrixCSC)
res = zeros(eltype(sp),sp.n,sp.n)
col = 1
@inbounds for i=1:length(sp.nzval)
while sp.colptr[col+1]<=i
col+= 1
end
col2 = 1
@inbounds for j in 1:length(sp.nzval)
if sp.rowval[j] != sp.rowval[i]
continue
end
while j>=sp.colptr[col2+1]
col2+= 1
end
res[col,col2] += sp.nzval[i]*sp.nzval[j]
end
end
res
end
On Wednesday, November 18, 2015 at 4:18:37 AM UTC+2, Dan wrote:
>
> It works, but it is slow (compared to `full(A'*A)`)
>
> On Wednesday, November 18, 2015 at 3:56:31 AM UTC+2, Dan wrote:
>>
>> And finally, getting rid of debug print.
>>
>> function AtrA(sp::SparseMatrixCSC)
>> res = zeros(eltype(sp),sp.n,sp.n)
>> col = 1
>> @inbounds for i=1:length(sp.nzval)
>> while sp.colptr[col+1]<i
>> col+= 1
>> end
>> col2 = 1
>> @inbounds for j in 1:length(sp.nzval)
>> if sp.rowval[j] != sp.rowval[i]
>> continue
>> end
>> while j>=sp.colptr[col2+1]
>> col2+= 1
>> end
>> res[col,col2] += sp.nzval[i]*sp.nzval[j]
>> end
>> end
>> res
>> end
>>
>>
>>
>>
>> On Wednesday, November 18, 2015 at 3:54:56 AM UTC+2, Dan wrote:
>>>
>>> Oops... a bug (which showed itself when matrix was sparser than my
>>> previous tests):
>>> [when columns where completely empty, the first `if` should have been a
>>> `while`]
>>>
>>> function AtrA(sp::SparseMatrixCSC)
>>> res = zeros(eltype(sp),sp.n,sp.n)
>>> col = 1
>>> @inbounds for i=1:length(sp.nzval)
>>> while sp.colptr[col+1]<i
>>> col+= 1
>>> end
>>> col2 = 1
>>> @inbounds for j in 1:length(sp.nzval)
>>> if sp.rowval[j] != sp.rowval[i]
>>> continue
>>> end
>>> while j>=sp.colptr[col2+1]
>>> col2+= 1
>>> end
>>> @show col,col2,sp.nzval[i],sp.nzval[j]
>>> res[col,col2] += sp.nzval[i]*sp.nzval[j]
>>> end
>>> end
>>> res
>>> end
>>>
>>>
>>>
>>>
>>> On Wednesday, November 18, 2015 at 3:35:30 AM UTC+2, Dan wrote:
>>>>
>>>> This is an O(nnz(A)^2) implementation. But, benchmarks will determine
>>>> the best solution.
>>>>
>>>> function AtrA(sp::SparseMatrixCSC)
>>>> res = zeros(eltype(sp),sp.n,sp.n)
>>>> col = 1
>>>> @inbounds for i=1:length(sp.nzval)
>>>> if sp.colptr[col+1]==i
>>>> col+= 1
>>>> end
>>>> col2 = 1
>>>> @inbounds for j in 1:length(sp.nzval)
>>>> if sp.rowval[j] != sp.rowval[i]
>>>> continue
>>>> end
>>>> while j>=sp.colptr[col2+1]
>>>> col2+= 1
>>>> end
>>>> res[col,col2] += sp.nzval[i]*sp.nzval[j]
>>>> end
>>>> end
>>>> res
>>>> end
>>>>
>>>>
>>>>
>>>>
>>>> On Wednesday, November 18, 2015 at 2:30:14 AM UTC+2, Dan wrote:
>>>>>
>>>>> Well, its not as simple as that because of everything being sparse
>>>>> (adding columns is not just a loop over sum). Actually, it would be
>>>>> interesting to see a clever solution.
>>>>>
>>>>> On Wednesday, November 18, 2015 at 2:14:39 AM UTC+2, Dan wrote:
>>>>>>
>>>>>> the columns of A'A when A is sparse are just sparse linear
>>>>>> combinations of columns of A. the coefficients being stored in the
>>>>>> column-based sparse matrix representation (internally a vector of
>>>>>> coefficients with pointers to new-column indices). so simply generating
>>>>>> column by column with a for loop should be as-quick-as-it-gets without
>>>>>> GPUing or some Fourier trickery.
>>>>>> did you already consider this? (i can try and make this description
>>>>>> into code)
>>>>>>
>>>>>> On Tuesday, November 17, 2015 at 10:26:21 PM UTC+2, Douglas Bates
>>>>>> wrote:
>>>>>>>
>>>>>>> More careful reading of the spmatmul function in
>>>>>>> base/sparse/linalg.jl provided a method directly downdates the dense
>>>>>>> square
>>>>>>> matrix C (i.e. does not form a sparse A'A then downdate) but I still
>>>>>>> need
>>>>>>> to form a transpose. If anyone knows of a fast way of forming A'A
>>>>>>> without
>>>>>>> creating the transpose I would be happy to hear of it.
>>>>>>>
>>>>>>
