Skip to content

Issue/Question about the output of lstsq #817

Closed
@loiseaujc

Description

@loiseaujc

Hej,

I'm really happy that the linalg module is on its way to being production-ready. I have already played a bit with the latest additions in v0.6.0 and observed something an undesirable (in my opinion) behaviour regarding the output of lstsq. Below is a MWE.

program main
  use iso_fortran_env, only: output_unit
  use stdlib_kinds, only: dp
  use stdlib_linalg, only: lstsq, solve_lstsq
  implicit none

  integer :: m, n
  real(dp), allocatable :: A(:, :), b(:), x(:), x_lstsq(:)

  ! Create the problem.
  m = 10 ; n = 5
  allocate(A(m, n)) ; call random_number(A)
  allocate(x(n)) ; call random_number(x)
  b = matmul(A, x)
 
  ! Solve the least-squares problem.
  x_lstsq = lstsq(A, b)

  write(output_unit, *) shape(A), shape(b), shape(x), shape(x_lstsq)
  write(output_unit, *) "Norm of the error :", norm2(x_lstsq(1:n)-x)

end program main

When running this code, the vector x_lstsq computed by lstsq is of size m rather than being of size n. I know that, by default, *gelsd takes A and b as input and overwrites the first n entries of b with the solution x. I believe this is quite misleading and can potentially cause issues if the user is directly using solve_lstsq. Calling solve_lstsq with x being an n-vector instead of an m-vector actually raises an error since the following test then fails

       if (lda<1 .or. n<1 .or. ldb<1 .or. ldb/=m .or. ldx/=m) then
          err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid sizes: a=',[lda,n], &
                                                           'b=',[ldb,nrhs],' x=',[ldx,nrhsx])
          call linalg_error_handling(err0,err)
          if (present(rank)) rank = 0
          return
       end if

in e.g. stdlib_linalg_s_solve_lstsq_one. Are there any reasons why we should keep the output vector as an m-vector rather than an n-vector ? From a mathematical and practical point of view, I would expect the following behaviour:

  • lstsq takes an m x n matrix A and an m-vector b as entry and returns an n-vector.
  • solve_lstsq takes as arguments m x n matrix A, an m-vector b and an n-vector x.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions