linearproj: Linear Projection onto a Subspace

In pracma: Practical Numerical Math Functions

Linear Projection onto a Subspace

Description

Computes the projection of points in the columns of B onto the linear subspace spaned by the columns of A, resp. the projection of a point onto an affine subspace and its distance.

Usage

  linearproj(A, B)

  affineproj(x0, C, b, unbound = TRUE, maxniter = 100)

Arguments

A	Matrix whose columns span a subspace of some R^n.
B	Matrix whose columns are to be projected.
x0	Point in R^n to be projected onto C x = b.
C, b	Matrix and vector, defining an affine subspace as C x = b
unbound	Logical; require all x >= 0 if unbound is false.
maxniter	Maximum number of iterations (if is unbound is false).

Details

linearproj projects points onto a linear subspace in R^n. The columns of A are assumed be the basis of a linear subspace, esp. they are required to be linearly independent. The columns of matrix B define points in R^n that will be projected onto A, and their resp. coefficients in terms of the basis in A are computed.

The columns of A need to be linearly independent; if not, generate an orthonormal basis of this subspace with orth(A). If you want to project points onto a subspace that is defined by A x = 0, then generate an orthonormal basis of the nullspace of A with null(A).

Technically, the orthogonal projection can be determined by a finite 'Fourier expansion' with coefficients calculated as scalar products, see the examples.

affineproj projects (single) points onto an affine subspace defined by A x = b and calculates the distance of x0 from this subspace. The calculation is based on the following formula:

p = (I - A' (A A')^{-1}) x0 + A' (A A')^{-1} b

Technically, if a is one solution of C x = b, then the projection onto C can be derived from the projection onto S = {C x = 0} with proj_C(x) = a + proj_S(x - a), see the examples.

In case the user requests the coordinates of the projected point to be positive, an iteration procedure is started where negative coordinates are set to zero in each iteration.

Value

The functions linearproj returns a list with components P and Q. The columns of P contain the coefficients – in the basis of A – of the corresponding projected points in B, and the columns of Q are the the coordinates of these points in the natural coordinate system of R^n.

affineproj returns a list with components proj, dist, and niter. proj is the projected point, dist the distance from the subspace (and niter the number of iterations if positivity of the coordinates was requested.).

Note

Some timings show that these implementations are to a certain extent competitive with direct applications of quadprog.

Author(s)

Hans W. Borchers, partly based on code snippets by Ravi Varadhan.

References

G. Strang (2006). Linear Algebra and Its Applications. Fourth Edition, Cengage Learning, Boston, MA.

See Also

nullspace, orth

Examples

#-- Linear projection --------------------------------------------------

# Projection onto the line (1,1,1) in R^3
A <- matrix(c(1,1,1), 3, 1)
B <- matrix(c(1,0,0, 1,2,3, -1,0,1), 3, 3)
S <- linearproj(A, B)
## S$Q
##           [,1] [,2] [,3]
## [1,] 0.3333333    2    0
## [2,] 0.3333333    2    0
## [3,] 0.3333333    2    0

# Fourier expansion': sum(<x0, a_i> a_i /<a_i, a_i>), a_i = A[ ,i]
dot(c(1,2,3), A) * A / dot(A, A)    # A has only one column

#-- Affine projection --------------------------------------------------

# Projection onto the (hyper-)surface x+y+z = 1 in R^3
A <- t(A); b <- 1
x0 <- c(1,2,3)
affineproj(x0, A, b)            # (-2/3, 1/3, 4/3)

# Linear translation: Let S be the linear subspace and A the parallel
# affine subspace of A x = b, a the solution of the linear system, then
#   proj_A(x) = a + proj_S(x-a)
a <- qr.solve(A, b)
A0 <- nullspace(A)
xp <- c(a + linearproj(A0, x0 - a)$Q)
## [1] -0.6666667  0.3333333  1.3333333

#-- Projection with positivity ----------------------- 24 ms -- 1.3 s --
s <- affineproj(x0, A, b, unbound = FALSE)
zapsmall(s$proj)                 # [1] 0 0 1
## $x     : 0.000000e+00 3.833092e-17 1.000000e+00
## $niter : 35

#-- Extended Example ------------------------------------------ 80 ms --
## Not run: 
set.seed(65537)
n = 1000; m = 100                       # dimension, codimension
x0 <- rep(0, n)                         # project (0, ..., 0)
A <- matrix(runif(m*n), nrow = m)       # 100 x 1000
b <- rep(1, m)                          # A x = b, linear system
a <- qr.solve(A, b)                     # A a = b, LS solution
A0 <- nullspace(A)                      # 1000 x 900, base of <A>
xp <- a+drop(A0 %*% dot(x0-a, A0))      # projection
Norm(xp - x0)                           # [1] 0.06597077

## End(Not run)

#-- Solution with quadprog ------------------------------------ 40 ms --
# D <- diag(1, n)             # quadratic form
# A1 <- rbind(A, diag(1, n))  # A x = b and
# b1 <- c(b, rep(0, n))       #   x >= 0
# n <- nrow(A)
# sol = quadprog::solve.QP(D, x0, t(A1), b1, meq = n)
# xp <- sol$solution

#-- Solution with CVXR ---------------------------------------- 50 ms --
# library(CVXR)
# x = Variable(n)                             # n decision variables
# objective = Minimize(p_norm(x0 - x))        # min! || p0 - x ||
# constraint = list(A %*% x == b, x >= 0)     # A x = b, x >= 0
# problem = Problem(objective, constraint)
# solution = solve(problem)                   # Solver: ECOS
# solution$value                              # 
# xp <- solution$getValue(x)                  # 

pracma documentation built on Nov. 10, 2023, 1:14 a.m.