Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implemented sparse QR decomposition #14

Merged
merged 5 commits into from
Oct 9, 2019
Merged

Implemented sparse QR decomposition #14

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
63ccb70
Added sparse QR
ZimmerA Nov 30, 2018
f3ff977
Got rid of "SparseEntry" type
ZimmerA Nov 30, 2018
8397205
Got rid of "SparseVector" type
ZimmerA Nov 30, 2018
ea5fad4
Refactored
ZimmerA Nov 30, 2018
bb35457
Refactored
ZimmerA Nov 30, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Got rid of "SparseVector" type
  • Loading branch information
@ZimmerA
ZimmerA committed Nov 30, 2018
commit 8397205e6a5d57e55b2110ff866cb47a40d13fa9
119 changes: 53 additions & 66 deletions src/FSharp.Stats/Algebra/LinearAlgebraServiceManaged.fs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,10 @@ namespace FSharp.Stats.Algebra
//open Microsoft.FSharp.Math
//open Microsoft.FSharp.Math. Bindings.Internals.NativeUtilities
open FSharp.Stats

open DotNetExtensions
/// This module is for internal use only.
module LinearAlgebraManaged =
open System

//let a = 1.

Expand Down Expand Up @@ -192,24 +193,6 @@ module LinearAlgebraManaged =
A.[l,k] <- A.[l,k] - 2.0 * v.[l-i] * v'A.[k-i]
v // Return reflection vector.

type SparseVector() =
member val entries = new ResizeArray<(int*float)>() with get,set
member m.Front
with get() = m.entries.[0]
and set(value) = m.entries.[0] <- value
member m.Add (entry:(int*float)) =
m.entries.Add(entry)
member m.IsEmpty =
m.entries.Count <= 0
member m.Swap (other:SparseVector) =
let temp = new ResizeArray<(int*float)> (other.entries)
other.entries <- new ResizeArray<(int*float)>(m.entries)
m.entries <- new ResizeArray<(int*float)>(temp)
member m.PopFront =
let value = m.entries.[0]
m.entries.RemoveAt(0)
value

/// QR factorization function for dense matrices
let QRDense (A:matrix) =
let (n,m) = matrixDims A
Expand Down Expand Up @@ -242,11 +225,11 @@ module LinearAlgebraManaged =

// Copies the row of a Sparsematrix into a Sparsevector
let CopyRow(nnz:int, offset:int, A:SparseMatrix<float>) =
let sparseVector = SparseVector()
let sparseVector = ResizeArray<(int*float)>()
for i=nnz-1 downto 0 do
if(not (A.SparseValues.[i+offset] = 0.0)) then
let entry = (A.SparseColumnValues.[i+offset], A.SparseValues.[i+offset])
sparseVector.entries.Insert(0, entry)
sparseVector.Insert(0, entry)
sparseVector

// Returns "value" with the sign of "sign"
Expand All @@ -266,7 +249,7 @@ module LinearAlgebraManaged =
Copysign(1.0/c, s)

// Apply the givens rotation on 2 Sparsevectors
let applyGivens( x : SparseVector, y : SparseVector) =
let applyGivens( x : ResizeArray<(int*float)>, y : ResizeArray<(int*float)>) =
let a = snd x.Front
let b = snd y.Front

Expand Down Expand Up @@ -298,100 +281,104 @@ module LinearAlgebraManaged =
let mutable q = 0
p <- p + 1 // Skip the first entry since we already handled it above

while (p < x.entries.Count && q < y.entries.Count) do
if (fst x.entries.[p] = fst y.entries.[q]) then
let xNew = c * snd x.entries.[p] - s * snd y.entries.[q]
let yNew = s * snd x.entries.[p] + c * snd y.entries.[q]
while (p < x.Count && q < y.Count) do
if (fst x.[p] = fst y.[q]) then
let xNew = c * snd x.[p] - s * snd y.[q]
let yNew = s * snd x.[p] + c * snd y.[q]
if (not (xNew = 0.0)) then
x.entries.[p] <- (fst x.entries.[p], xNew)
x.[p] <- (fst x.[p], xNew)
p <- p + 1
else
x.entries.RemoveAt(p)
x.RemoveAt(p)
if ( not (yNew = 0.0)) then
y.entries.[q] <- (fst y.entries.[q], yNew)
y.[q] <- (fst y.[q], yNew)
q <- q + 1
else
y.entries.RemoveAt(q)
else if(fst x.entries.[p] < fst y.entries.[q]) then
let k = fst x.entries.[p]
let xNew = c * snd x.entries.[p]
let yNew = s * snd x.entries.[p]
x.entries.[p] <- (fst x.entries.[p] , xNew)
y.RemoveAt(q)
else if(fst x.[p] < fst y.[q]) then
let k = fst x.[p]
let xNew = c * snd x.[p]
let yNew = s * snd x.[p]
x.[p] <- (fst x.[p] , xNew)
p <-p + 1
y.entries.Insert(q,(k, yNew))
y.Insert(q,(k, yNew))
q <- q + 1
else
let k = fst y.entries.[q]
let xNew= -s * snd y.entries.[q]
let yNew= c * snd y.entries.[q]
x.entries.Insert(p,(k, xNew))
let k = fst y.[q]
let xNew= -s * snd y.[q]
let yNew= c * snd y.[q]
x.Insert(p,(k, xNew))
p <- p + 1
y.entries.[q] <- (fst y.entries.[q], yNew)
y.[q] <- (fst y.[q], yNew)
q <- q + 1
if(p < x.entries.Count) then
while(p < x.entries.Count) do
let k = fst x.entries.[p]
let xNew = c * snd x.entries.[p]
let yNew = s * snd x.entries.[p]
x.entries.[p] <- (fst x.entries.[p], xNew)
if(p < x.Count) then
while(p < x.Count) do
let k = fst x.[p]
let xNew = c * snd x.[p]
let yNew = s * snd x.[p]
x.[p] <- (fst x.[p], xNew)
p <- p + 1
y.entries.Insert(q,(k, yNew))
y.Insert(q,(k, yNew))
q <- q + 1
else if(q < y.entries.Count) then
while(q < y.entries.Count) do
let k = fst y.entries.[q]
let xNew = -s * snd y.entries.[q]
let yNew = c * snd y.entries.[q]
x.entries.Insert(p,(k, xNew))
else if(q < y.Count) then
while(q < y.Count) do
let k = fst y.[q]
let xNew = -s * snd y.[q]
let yNew = c * snd y.[q]
x.Insert(p,(k, xNew))
p <- p + 1
y.entries.[q] <- (fst y.entries.[q], yNew)
y.[q] <- (fst y.[q], yNew)
q <- q + 1

encode(c, s)

// Run Sparse QR
let m = A.NumRows
let n = A.NumCols
let mutable Q = new ResizeArray<SparseVector>()
let mutable R = new ResizeArray<SparseVector>()
let mutable Q = new ResizeArray<ResizeArray<(int*float)>>()
let mutable R = new ResizeArray<ResizeArray<(int*float)>>()

for i = 0 to m - 1 do
Q.Add(SparseVector())
R.Add(SparseVector())
Q.Add(ResizeArray<(int*float)>())
R.Add(ResizeArray<(int*float)>())

for a = 0 to m-1 do
let row = CopyRow(A.SparseRowOffsets.[a+1] - A.SparseRowOffsets.[a], A.SparseRowOffsets.[a], A)
let mutable row = CopyRow(A.SparseRowOffsets.[a+1] - A.SparseRowOffsets.[a], A.SparseRowOffsets.[a], A)
let mutable q = 0
while((not row.IsEmpty) && fst row.Front < a && fst row.Front < n) do
let j = fst row.Front

if(R.[j].IsEmpty || fst R.[j].Front > j) then
R.[j].Swap(row)
Q.[a].entries.Insert(q, (j, 1.0))
let temp = new ResizeArray<(int*float)> (row)
row <- new ResizeArray<(int*float)>(R.[a])
R.[a] <- new ResizeArray<(int*float)>(temp)
Q.[a].Insert(q, (j, 1.0))
q <- q + 1
()
else
let ret = applyGivens(R.[j], row)
Q.[a].entries.Insert(q, (j, ret))
Q.[a].Insert(q, (j, ret))
q <- q + 1
()
if (a < n) then
R.[a].Swap(row)
let temp = new ResizeArray<(int*float)> (row)
row <- new ResizeArray<(int*float)>(R.[a])
R.[a] <- new ResizeArray<(int*float)>(temp)
()

let mutable x = -1

let RSeq = seq{ for row in R do
x <- x+1
for entry in row.entries do
for entry in row do
yield(x, fst entry, snd entry)}

let ROut = Matrix.initSparse m n RSeq

x <- -1
let QSeq = seq{ for row in Q do
x <- x+1
for entry in row.entries do
for entry in row do
yield(x, fst entry, snd entry)}

let QOut = Matrix.initSparse m n QSeq
Expand Down
15 changes: 15 additions & 0 deletions src/FSharp.Stats/DotNetExtensions.fs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,15 @@
module DotNetExtensions

open System.Collections.Generic

[<Extension>]
type List<'T> with
member m.Front
with get() = m.[0]
and set(value : 'T) = m.[0] <- value
member m.PopFront =
let value = m.[0]
m.RemoveAt(0)
value
member m.IsEmpty =
m.Count <= 0
1 change: 1 addition & 0 deletions src/FSharp.Stats/FSharp.Stats.fsproj
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
<EmbeddedResource Remove="Scripts\**" />
</ItemGroup>
<ItemGroup>
<Compile Include="DotNetExtensions.fs" />
<Compile Include="AssemblyInfo.fs" />
<Compile Include="Ops.fs" />
<Compile Include="Random.fs" />
Expand Down