Randomized Linear Algebra

This notebook aims to show the utility of the Lighton OPU for randomized linear algebra tasks. See https://arxiv.org/abs/2104.14429 for more details.

import lightonml
from lightonml import OPU
import numpy as np
from lightonml.encoding import base
from tqdm.notebook import tqdm

opu = OPU(n_components = 2000)

The lighton OPU internally does random projections of binary vectors. However, for randomized linear algebra we generally need linear projections of floating point inputs.

To make a linear proections with binary inputs, it is sufficient to call opu.linear_transform.

To handle floating point inputs, we use the encoders from lightonml.encoding. In particular we will use SeparatedBitPlanEncoder which will split the input into mutliple bit planes and then transform them one after the other. SeparatedBitPlanEncoder will then recombine it, preserving the linear projection.

encoder = base.SeparatedBitPlanEncoder
decoder = base.SeparatedBitPlanDecoder

We will now showcase three applications of randomized linear algebra Trace Estimation, Triangle Estimation and Randomized SVD.

For trace estimation we consider Hutchinson’s estimator \(Tr(A) \approx Tr(R^\top AR)\)

Where \(R\) is an appropriately scaled normal random matrix.

We can achieve this with the following code:

def trace_estimation(A, m):
    n = A.shape[0]
    opu.fit1d(A)
    s = np.random.normal(0,1,(100, n))

    opu_var = opu.linear_transform(s, encoder, decoder).var()
    AR   = opu.linear_transform(A,    encoder, decoder)[:, :m]
    RTAR = opu.linear_transform(AR.T, encoder, decoder)[:, :m]
    return RTAR.trace() / m * n / opu_var

For triangle estimation we combine Hutchinson’s estimator with sketched matrix multiplication, which gives the following expression:

\(Tr(A^3) \approx Tr(R^\top AR)^3\)

def triangle_estimation(A, m):
    n = A.shape[0]
    opu.fit1d(A)
    s = np.random.normal(0,1,(100, n))

    opu_var = opu.linear_transform(s, encoder, decoder).var()
    AR   = opu.linear_transform(A,    encoder, decoder)[:, :m]
    RTAR = opu.linear_transform(AR.T, encoder, decoder)[:, :m]
    return np.linalg.matrix_power(RTAR/ opu_var*n/m,3).trace()

Finally, following the procedure in Halko et al, the singular value decompositon can also be approximated using RandNLA. This procedure can also be accelearated using the OPU.

def randomized_svd(A, m, num_iter = 0):
    Q = opu.linear_transform(A, encoder, decoder)[:, :m]
    opu.fit1d(A)
    Q, _ = np.linalg.qr(Q)
    for it in range(num_iter): # power_iteration
        Q = A@Q
        Q, _ = np.linalg.qr(Q)
        Q = (Q.T@ A).T
        Q, _ = np.linalg.qr(Q)
    U, s, Ra = np.linalg.svd(A@Q, full_matrices=False)
    Va = Ra@Q.T
    return U[:, :m], s[:m], Va[:m, :]

def relative_error(approximtion, correct):
    return np.linalg.norm(approximtion-correct)/ np.linalg.norm(correct)

Testing the functions

Trace Estimation

A = np.random.uniform(0,1,(2000,2000))
opu.fit1d(A)
m = 200
estimated_trace= trace_estimation(A, m)
real_trace = np.trace(A)

print("Estimated Trace: {:.2f}\nReal Trace:      {:.2f}\nRelativeError:   {:.2f}".format(estimated_trace, real_trace, relative_error(real_trace, estimated_trace)))

Estimated Trace: 1045.87
Real Trace:      990.05
RelativeError:   0.05

Triangle Estimation

import networkx as nx
graph = nx.generators.random_graphs.barabasi_albert_graph(2000, 500)
A = np.array(nx.linalg.graphmatrix.adjacency_matrix(graph).todense())
opu.fit1d(A)
m = 200
estimated_tris= triangle_estimation(A, m)
real_tris = np.trace(np.linalg.matrix_power(A,3))

print("Estimated Triangles: {:.2f}\nReal Triangles:      {:.2f}\nRelativeError:       {:.2f}".format(estimated_tris, real_tris, relative_error(estimated_tris, real_tris)))

Estimated Triangles: 783703186.68
Real Triangles:      586469286.00
RelativeError:       0.34

Randomized SVD

A = np.random.normal(0,1,(2000,2000))
opu.fit1d(A)
m = 1000
# randomized svd.
Ur, Sr, Vr = randomized_svd(A, m, num_iter = 2)
# exact (truncated) svd.
Ue, Se, Ve = np.linalg.svd(A)

relative_error_rsvd = relative_error(Ur@np.diag(Sr)@Vr, A)
relative_error_tsvd = relative_error(Ue[:, :m]@np.diag(Se[:m])@Ve[:m,:], A)

print("Truncated SVD:  {:.3f}\nRandomized SVD: {:.3f}".format(relative_error_tsvd, relative_error_rsvd))

Truncated SVD:  0.326
Randomized SVD: 0.351